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A majority of methods from dynamical systems analysis, especially those in applied settings, rely on 
Poincare's geometric picture that focuses on "dynamics of states" . While this picture has fueled our 
field for a century, it has shown difficulties in handling high-dimensional, ill-described, and uncertain 
systems, which are more and more common in engineered systems design and analysis of "big data" 
measurements. This overview article presents an alternative framework for dynamical systems, based 
on the "dynamics of observables" picture. The central object is the Koopman operator: an infinite- 
dimensional, linear operator that is nonetheless capable of capturing the full nonlinear dynamics. 
The ffist goal of this paper is to make it clear how methods that appeared in different papers and 
contexts all relate to each other through spectral properties of the Koopman operator. The second 
goal is to present these methods in a concise manner in an effort to make the framework accessible 
to researchers who would like to apply them, but also, expand and improve them. Finally, we aim 
to provide a road map through the literature where each of the topics was described in detail. We 
describe three main concepts: Koopman mode analysis, Koopman eigenquotients, and continuous 
indicators of ergodicity. For each concept we provide a summary of theoretical concepts required 
to define and study them, numerical methods that have been developed for their analysis, and, 
when possible, applications that made use of them. The Koopman framework is showing potential 
for crossing over from academic and theoretical use to industrial practice. Therefore, the paper 
highlights its strengths, in applied and numerical contexts. Additionally, we point out areas where 
an additional research push is needed before the approach is adopted as an off-the-shelf framework 
for analysis and design. 
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A majority of methods from dynamical sys- 
tems analysis, especially those in applied set- 
tings, rely on Poincare's geometric picture 
that focuses on "dynamics of states". While 
this picture has fueled our field for a century, 
it has shown difficulties in handling high- 
dimensional, ill-described, and uncertain sys- 
tems, which are more and more common in 
engineered systems design and analysis of 
"big data" measurements. This overview ar- 
ticle presents an alternative framework for 
dynamical systems, based on the "dynam- 
ics of observables" picture. We present an 
overview of several approaches to studying 
dynamical systems using the Koopman op- 
erator, which holds promise to resolve these 
issues. The dynamics are analyzed by looking 
at evolutions of functions on the state space, 
rather than directly at state space trajecto- 
ries. The evolution can be understood by 
expanding the function into a basis of eigen- 
functions of the Koopman operator. The first 
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approach is based on the Koopman modes, 
which generalize linear mode analysis from 
linear systems to nonlinear systems, while 
preserving global nonlinear features of the 
system, unlike, e.g., linearizations based on 
Taylor- and Fourier- expansions. The sec- 
ond approach identifies coherent structures 
in fiows. An equivalence relation between 
points in the state space can be defined us- 
ing spectral properties of the Koopman oper- 
ator, where equivalent points correspond to 
initial conditions that behave statistically the 
same with respect to any observable. The 
third approach we present introduces contin- 
uous quantifications of ergodicity and mix- 
ing, concepts existing in ergodic theory that 
are traditionally treated as binary notions. 
Throughout the paper, we highlight exam- 
ples from the literature using each of these 
concepts. Examples are taken from diverse 
areas such as fiuid mechanics, fiuid mixing, 
energy efficiency of buildings, power systems, 
and Unmanned Aerial Vehicle path-planning 
for search-and-rescue. A common trait of all 
the methods is that they do not require ac- 
cess to an analytical model of the system; the 
spectral properties of the Koopman operator 
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can be constructed from measured or simu- 
lated data. 



I. INTRODUCTION 

Currently, dynamical systems analysis and design 
primarily uses the geometric picture, as put forth 
by Poincare in his work on the three body prob- 
lem. Much of the framework is built around notions 
from differential geometry, trajectories and invariant 
manifolds. Such an approach has met with success 
in a variety of settings and, at this point, one hardly 
needs to justify the use of geometric theory when 
working on a particular problem. 

However, the geometric viewpoint is ill-suited to 
many of the situations that are of interest in real 
systems. For example, for systems possessing hy- 
perbolic regimes, the unstable manifolds give rise 
to locally exponentially divergent trajectories. Any 
noise or uncertainty in the system will lead to multi- 
ple possible trajectories for an initial condition, with 
the width of the set trajectories initially expanding 
exponentially. In such cases, questions about the be- 
havior of a specific trajectory are difficult to answer. 

Systems with a large number of dimensions can 
be problematic as well, since many of the geometric 
arguments are only valid in a low number of dimen- 
sions, e.g., Bendixson's Criterion for determining the 
non-existence of periodic orbits in the plane. In some 
cases, these arguments can be extended, with diffi- 
culty, to an arbitrary number of dimensions. Even 
in these cases, however, a practical implementation 
limits them to a moderate number of dimensions. 
To handle high-dimensional systems, special sym- 
metries or other conditions are required in order 
to effectively reduce the dimension to a manage- 
able size. Furthermore, without access to explicit 
ODEs, even basic geometric analysis is difficult to 
apply. If dynamical systems theory is to become 
an important field in the context of pressing prob- 
lems such as "big data" , tools need to be developed 
that are capable of handling high-dimensional, un- 
certain, and ill-described systems, as well as systems 
for which past time-evolution data is available, but 
for which no simple mathematical description can be 
determined.'^ 

This article presents a viewpoint that is at the 
intersection of applied ergodic theory and opera- 
tor theory. These two fields can be used in ap- 
plied settings to analyze and design dynamical sys- 
tems, with many of the aforementioned difficulties 
being handled with a certain amount of elegance. 
In fact, when we study dynamical systems through 
certain linear operators, the full nonlinear dynam- 



ics can be captured within a linear setting. This 
linear setting allows the power of spectral analy- 
sis to be brought to bear on a (nonlinear) problem 
without sacrificing any information as required by 
other linearization techniques. Contrast this with 
the traditional spectral approach that only deter- 
mines geometry locally in the state space. Addi- 
tionally, in theory, the operator-theoretic approach 
works equally well whether the original state-space is 
low- or high-dimensional; the same techniques apply 
to both cases. The framework is also well suited to 
studying noisy systems because the primary object 
of interest is no longer the trajectory. Finally, and 
perhaps most significantly, the operators involved 
can be constructed, approximated, or analyzed using 
only simulation or experimental measurement data. 
This allows a certain black-box approach to the anal- 
ysis which is quite useful in real problems where the 
practitioner may not have full knowledge of the sys- 
tem's internals. 

As with any technique, however, there is a tradeoff 
in order to gain the above advantages. The operator- 
theoretic picture has no immediate connection to our 
physical intuition, making its meaning more diffi- 
cult to comprehend. One's viewpoint must change 
from considering the evolution of points in the state 
space to considering the evolution of functions. Ad- 
ditionally, the new approach is inherently infinite- 
dimensional, even when the state space is finite- 
dimensional. This sacrifice is what allows the full 
information of a nonlinear system to be contained 
within a linear setting. Because of this, the imple- 
mentation of any approximation is a more delicate 
issue. Finally, the associated numerical techniques 
are underdeveloped. Most of our approaches employ 
direct computations which are little more than nu- 
merical implementations of proofs. 

The two main candidates for the study of systems 
via operators are the Koopman operator and the 
Perron-Frobenius operator. In appropriate function 
spaces, they are duals to each other, so theoreti- 
cally, there should not be any distinction in working 
with one as opposed to the other. However, as men- 
tioned previously, we must always include applied 
considerations. Questions arise such as how do we 
construct or represent the chosen operator from the 
problem description and given data? How well does 
a finite approximation represent the ideal theoreti- 
cal picture? What part of intuition gained is due to 
numerical artifacts and what is real? 

The Perron-Frobenius operator represents a "dy- 
namics of densities" picture; it looks at groups of 
trajectories. One can think of this as watching the 
evolution of a mass distribution under the action of 
a fiow. From a numerical perspective, construction 
of the operator relies on selecting a set of initial con- 
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ditions and simulating forward for only a short time 
period, thus avoiding the compounding of numerical 
time-integration errors. Due to these short bursts, 
transient dynamics can be captured very well. How- 
ever, much attention has been focused on computing 
invariant densities,'^ which are infinite-time objects, 
through approximating the Perron-Frobenius oper- 
ator by a Markov chain. The number of simulated 
initial conditions is dictated by the need to sam- 
ple the region of interest well. In high-dimensions, 
both short- and long-time dynamics simulations re- 
quire a mesh on the entire space. This can be true 
even in the case of a low-dimensional attractor. If 
we have a priori knowledge of the low-dimensional 
subspace the attractor lives in, then the mesh size 
can be restricted. However, for an arbitrary system, 
this knowledge may not be initially available, thus 
requiring the full mesh. 

On the other hand, the Koopman operator 
presents a picture for the "dynamics of observables" . 
The difference in viewpoints between the Perron- 
Frobenius and Koopman operators is similar to the 
Eulerian versus the Lagrangian viewpoint in fiuid 
mechanics, with the Koopman picture correspond- 
ing to the Lagrangian viewpoint. What is meant is 
that measurements are made along trajectories. For 
the Koopman operator, the numerical construction 
relies on potentially fewer initial conditions, but re- 
quires longer run-times, which is more suitable to 
physical experiments. For example, when testing 
a jet engine, it is started from a relatively small 
number of initial conditions and run it over a long 
time rather than preparing thousands of initial con- 
ditions and running the engine for a few seconds for 
each initial condition. Due to the long run-times re- 
quired, the asymptotics are well-understood. How- 
ever, more research is needed to understand the tran- 
sients. 

To visualize high-dimensional dynamical systems, 
we often restrict our attention to one, or a few, two- 
dimensional cross-sections in the state space and 
look at the invariant structures intersecting that 
slice. With the Perron-Frobenius operator, it is dif- 
ficult to directly compute invariant densities on the 
slice of interest, since, in principle, it requires a com- 
putation of the invariant density for the entire state 
space as an intermediate step. For the Koopman op- 
erator, invariant objects are attached to initial con- 
ditions, making it well-suited to visualizing struc- 
tures on an arbitrary 2D cross-section in the state 
space. Initial conditions can be easily prepared on 
the slice and the invariants directly computed. In 
such cases, the number of initial conditions required 
to understand the dynamics is significantly reduced. 

While operator methods, and specifically the 
Koopman and Perron-Frobenius operators, have 



much potential to deal with applied problems, these 
methods are all but absent from the a pplied a nd 
industrial settings, with due except ions .1^^^^^^ To 
speed up the adoption of operator techniques in 
these domains, any new methodology needs to be 
able to leverage already existing data, instead of 
proposing both a new methodology and a new way 
to collect data. The "dynamics of observables" per- 
spective, and specifically the Koopman operator, is 
used as it deals with measurements, i.e., observables, 
which are well-understood both theoretically and 
computationally. On the other hand, the Perron- 
Frobenius techniques would require working with 
representations of densities, which are often singu- 
lar, especially in well-behaved engineered systems. 

In this paper, we intend to describe three con- 
cepts, all under the umbrella of "dynamics of ob- 
servables", that show how this theory can be made 
useful for analysis and design. Contributions can be 
split between theoretical and applied contributions. 



Theoretical contributions 

1. Dynamical evolution of a system can be stud- 
ied by looking at what is termed Koopman mode 
analysis. The concept is similar to normal mode 
analysis familiar from linear vibration theory. Koop- 
man mode analysis starts with a choice of a set of 
linearly independent observables, or equivalently a 
vector-valued observable. The Koopman operator 
U is then analyzed through its action on the sub- 
space spanned by the chosen observables. The ob- 
servables are decomposed into projections onto the 
eigenspaces of [/, and the evolution is a sum of terms 
composed of a product of three terms: i) a part that 
is time-dependent and is determined by the eigen- 
value (or frequency) associated with the eigenspace; 
ii) an eigenfunction of /7, which is a function of 
the initial conditions; iii) the vector of the coeffi- 
cients of the projection of the observables onto the 
eigenspaces, with the coefficients only being func- 
tions of the chosen observables. In this way, spec- 
tral analysis can be performed on nonlinear systems. 
This analysis is also used for model reduction. ^ 

2. The notions of the ergodic quotients and eigen- 
quotients allow the Koopman operator to be used 
for the extraction and analysis of invariant and pe- 
riodic structures in the state space.^ The points in 
the state space are grouped into invariant sets using 
level sets of eigenfunctions of the Koopman oper- 
ator. Instead of set-theoretic framework, this ap- 
proach is lifted to the analytical setting using the 
ergodic quotient and eigenquotient formalism. The 
eigenquotients are studied as subsets of particular 
Sobolev spaces, where their geometry gives insight 
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into the structure of the state space, in spirit similar 
to analysis of Hamiltonian systems via Morse theory 
of the associated energy functions. 

3. In the standard interpretation of ergodic the- 
ory, mixing and ergodicity are treated as binary con- 
cepts: a system is either mixing/ergodic or it is not. 
Both mixing and ergodicity can be formulated using 
spectral invariants of the Koopman operator. From 
a finite-time evolution of an arbitrary system, we 
can quantify how close its spectral invariants are 
to the "ideal" case, e.g., mixing or ergodic, and in 
this way formulate continuous indicators of ergod- 
icity and mixing. The ergodicity defect^^ and the 
mixing norm^^ are examples of continuous indicators 
that extend the corresponding binary notions. Such 
a relaxation brings the concepts of ergodicity and 
mixing into an engineering context, allowing, e.g., 
the use of the indicators as optimization criteria. We 
present a unified explanation of the co ncepts t hat 
have previously appeared in liter at ure.'^^M^^'^^EHl 



Numerical techniques and applications 

Numerical computation of the objects in the the- 
ory uses elements from three different areas. Fourier 
analysis based methods are useful for computing 
Koopman modes for dynamics on the attractor in 
addition to being essential for the construction of the 
eigenquotient. A variant of the standard Arnoldi al- 
gorithm based on companion matrices is also useful 
for computing part of the spectrum of the Koopman 
operator, a basic element of Koopman mode anal- 
ysis. This variant does not require an explicit rep- 
resentation of the operator and only requires data, 
sequences of vectors coming from either simulations 
or experiments. 

Koopman mode analysis has seen applications in 
flui ds mechan ics to extract spatial structures for the 
Aq-^ [81M1 47 | 49 | Koopman modes have also found appli- 
cations in the analysis of coherency and instabilities 
for power systems^^ and in the field of building en- 
ergy efficiency where they h ave be en used for model 
validation and data analysis .^^^^ 

To compute the eigenquotients numerically, a 
set of observables is averaged along trajectories 
started at different initial conditions, obtaining a 
finite-dimensional representation of any eigenquo- 
tient. Such representations are a nalyz ed with the 
aid of a diffusion maps algorithm 5^^^ which com- 
putes a change of coordinates, the diffusion modes, 
for the eigenquotient. In the limit where infinitely 
many initial conditions were simulated for infinite 
time, such coordinate change would render any con- 
sequent analysis independent of the choice of the 
observables averaged during the computation. The 



scale-ordering of diffusion coordinates makes it prac- 
tical to obtain a low-dimensional approximation of 
the eigenquotients. 

The averaging along trajectories is used again to 
formulate continuous indicators for ergodicity and 
mixing, where the rate of approach of the averages 
along finite-time trajectories to the infinite limit is 
indicative of the underlying dynamics. Based on 
such indicators, the dynamics of the fiows can be 
designed to match a particular statistical behavior, 
with applications in path-planning for vehicles^^ and 
mixing of fiuids on micro-scales.!^ 

II. NOTATION AND TERMINOLOGY 

We start off by fixing some notation and terminol- 
ogy. Let us denote the state space by M and define 
dynamics on it by the iterated map T : M ^ M. 
Note that the set M can be an arbitrary set (pos- 
sessing no structure) and T can be an arbitrary map 
on this set. Then the abstract dynamical system is 
specified by the couple (M, T). Note that standard 
texts on ergodic theory study a specific case when 
M is a measurable space, with a a- algebra 05, and 
T is 53-measurable. Additionally, transformation T 
is typically assumed to be measure-preserving, i.e., 
there exists a measure /i, the invariant measure, such 
that for any 5 G ^ 

^i{S) = ^Ji{T-^S), (1) 

with T~^S understood as the pre-image of S. The 
measure fi does not necessarily have a density func- 
tion associated with it. For the formulation of the 
theory in this paper, we do not require the measur- 
able framework, although when answering more spe- 
cific questions, we might restrict ourselves to it, as 
it is the one most commonly encountered in applied 
dynamical systems. 

We will be concerned with the behavior of observ- 
ables on the state space. To this end, we define an 
observable to be a function / : M ^ C, where / 
is an element of some function space For now, 
it is not necessary to specify any structure for 
A concrete interpretation of an observable is that 
of a sensor probe for the dynamical system in ques- 
tion; we can access information about the system via 
the evolution of the observable's values. Instead of 
tracking the trajectory {p, T(p), T^(p), . . . }, we now 
track the trace {/(p), f {T {p)) J {T\p)) , . . . }. The 
description of the dynamics can then be concisely 
written down in the form of state and output equa- 
tions, familiar to control theorists: 

Vn f{Pn)- 
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The dynamical systems community mainly focuses 
on state space trajectories {pn}, while the control 
systems community usually studies systems that 
have an additional input or disturbance terms, with 
the functions T and / taking particular forms that 
are common in engineered systems. 

We define the (discrete-time) Koopman operator, 
t/T : ^ ^ as 



[f/T/](p)=/(T(p)), 



(3) 



i.e., it is a composition, Ut/ = / o T, of the observ- 
able / and the iterated map T. When it is obvious 
which transformation gives rise to the Koopman op- 
erator, we will drop the dependence on T from the 
notation and write U instead of Ut- When is a 
vector space, is a linear operator. 

When M is a finite set, is a finite-dimensional 
operator and can be represented by a matrix. How- 
ever, when M is finite- or infinite-dimensional, U is 
generally infinite dimensional. Much of the time, we 
only have access to a particular collection of observ- 
ables {/i, . . . , fx} C J-*; these could be physically 
relevant observables arising naturally from the prob- 
lem, or a (subset of a) function basis for T. We can 
extend the Koopman operator to this larger space in 
the natural way: If F = (/i, . . . , /k)^ ^ then 
Uk ' J^^ is defined as 



[UkF]{p) :-- 



'[UhM 



(4) 



Hence Ur = (8)f" With an abuse of notation, we 
generally write Uk as U. The space J^^ is the space 
of C^- valued observables on M. In this context, 
is referred to as the output space. More generally, we 
can consider vector valued observables, F : M ^ V, 
where V is some vector space. For example, when 
analyzing the heat equation on a periodic box B, the 
state space can be regarded as the sequence space 
of Fourier coefficients and an observable F : M ^ 
can be regarded as mapping between a 
sequence of Fourier coefficients (the state space M) 
and a temperature distribution on B (the real- valued 
space L^(B,(ix)). We will revisit this setup in more 
detail in Example [4j 

The above notion of the Koopman operator was 
defined in the context of discrete-time dynamical 
systems. Often though, working in the discrete-time 
setting poses an unneeded restriction; in many sys- 
tems, the natural formulation of the dynamics is 
with respect to a continuous time variable. The 
Koopman operator can be extended to deal with 
continuous- time dynamical systems, or even more 
generally, event-based dynamical systems. 



Assume we have the continuous-time dynamical 
system p = T{p). In this context, there is not just 
the Koopman operator, but a semigroup of operators 
{U*}teR+ given by a generator U. We call the semi- 
group {U^} the Koopman semigroup. We explicitly 
define the action of the semigroup on the observable 
/ G as 



[U'f]{p) = f{'^'{p)). 



(5) 



Here ^^{p) = ^{p^t) is the ffow map that takes an 
initial condition p G M and maps it to the solution 
at time t of the initial value problem (IVP) having 
initial condition p{0) = p; i.e., for a fixed po G M, 
the trajectory {^{po^ ')}t>o is a solution of the IVP 
p = T{p)^ p{0) = pq. The generator of the Koopman 
semigroup is defined by 



[Uf] := lim 



u'f-f 

t ' 



(6) 



where the limit is taken in the strong sense.^ 

The following examples describe the above con- 
cepts in certain simple, concrete cases. 

Example 1 (Cyclic group). Let M = {e, a, a^} be a 
cyclic group of order 3 (a^ = e). Define T : M ^ M 
by T{p) = a - p. Hence the entire state space is a 
periodic orbit of period 3. Let be the C- valued 
functions on M. Clearly, the space of observables 
is C^. Let /i,/2,/3 be the indicator functions on 
e^a^o?^ respectively: 



hip) 



/2(P) 



1, 


p = 


e, 


0, 


p^ 


e 


1, 


p = 


a, 


0, 


p^ 


a 


1, 


p = 


a" 


0, 


p^ 


a? 



(7) 



2 • 



These form a basis for The action of the Koop- 
man operator on this basis is given as 



[ufiW) 

[Uf2]{p) 

[Uh]{jp) 



fi{a-p) 
f2{a-p) 
h{a-p) 



h{p). 
flip), 
/2(P). 



(8) 



For an arbitrary observable f ^ T given by / = 
Ci/i -\- €2/2 + C3/3, with Q G C, we have that 

Uf = C1/3 + C2/1 +C3/2. 

Then, the matrix representation of U in the 
{/i, /2, /3} basis is given by 



U 



Cl 




"0 1 0" 




Cl 


C2 




1 




C2 


C3_ 




1 




_C3_ 



(9) 
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In this case, ([9| gives the fuh action of the Koopman 
operator on T . □ 

Example 2 (Linear, diagonahzable systems). Let 
M = and define T : M ^ M by 



(10) 



where x = (xi, . . . , x^^)""" G M and /i^ G M. Let 
be the space of C- valued functions on W^. Let 
{6i,...,6c^} C M be a basis for M and define 
fi{x) = (bi^x)^ where (•, •) is the inner product on 
R^. The action of U : J' ^ J' on fi is 



[[//,](x) = (6„T(x)) = [6,,i,. 



Ml 

/i2 





(11) 

Let J^"^ = ®{j^ and define [/^ on J"^ as in Q. Then 
for F = (/i,...,/,)T, 



[UdF]{x) 



hdA 



Ml 

112 





Md 



Xd 



(12) 



Note that ( |T2| ) is just the action of the Koop- 
man operator on the particular observable F = 
(/i, • • • , /d)^, not the full action of the Koopman op- 
erator on the entire observable space, J^. The main 
point is that the Koopman operator is reducible, 
provided T leaves some subspace of T invariant. 

As a special case, we can take the 6i's to be the 
canonical basis vectors having zeros everywhere, ex- 
cept for a 1 in the i^^ entry. In this case, the func- 
tions {fi} represent the canonical projections onto 
the coordinates of x; i.e. fi{x) = Xi. Then the ac- 
tion of the Koopman operator with respect to these 
particular observables is given as 



[UdF]{x) 



Ml 
M2 





Md 



Xi 



Xd 



(13) 



□ 



The next example shows that the Koopman op- 
erator formalism can easily handle state spaces that 
are mixtures of discrete and continuous domains. 



Example 3 (Mixed state space). Let = [0, 2it) x 
[0, 27r) and G = {0, 1, 2}. Define M = T x G and let 
the dynamics be given by 



^/c+l 

<5fe+i 



I, 



K sin Ok, mod 27r 

4+1, mod 27T (1"^) 
1, mod 3 



where pj. = {h,Ok) ^ T^, Sk e G, and K > 0. 
The dynamics are given by the standard map cycling 
between the unperturbed, shear-fiow case and two 
perturbed cases. The "perturbation dynamics" are 
driven by a group action. 

Let Tt2 = {/ : ^ C} be the set of all func- 
tions mapping the torus into the complex numbers. 
We do not assume that the functions in J^j2 have 
any type of regularity or algebraic properties; for the 
moment they are completely arbitrary. Similarly, let 
J~G = : G ^ C} be the set of all functions from 
the group G into the complex numbers. One possi- 
ble choice for the space of observables on M is the 
set T = {h = g ■ f \ f G Tj2^g G J^g}- Hence ob- 
servables on M are pointwise products of functions 
on and G and map the mixed state space into C. 
The Koopman operator can easily be defined as 

[Uh]{pk, Sk) = g{sk+i) • /(P/e+l). 

Another possible choice for T could be the set of all 
the observables that are functions of only / and 0. 
This a a subset of the previous choice by taking g 
to be a constant function. This is a natural choice 
in the case that the "perturbation dynamics" (the 
dynamics on G), cannot be measured. 

□ 

Example 4 (Partial differential equations). Con- 
sider the 2D heat equation on B = [— ^, |] x [— ^, ^] 
with periodic boundary conditions: 



(15) 



Assuming u^V^u G I/^(B,(ix), u{x^t) can be ex- 
panded in a trigonometric basis: 



uix,t)= ^ a,-(t)e*2-J'-^ 



(16) 



where j x is the dot product of j and x. A Galerkin 
projection onto this basis yields 



aj{t) = -A7r^c^\\j\\laj{t). 



(17) 



Thus, we have the continuous-time, infinite- 
dimensional analogue of example |2j We could pro- 
ceed with exhibiting the Koopman semigroup or in- 
duce a discrete-time evolution from the continuous- 
time flow. In the latter case, fix a time step h > 



7 



and get 



ajit„+i) = exp{-4Tr c \\j\\2h) aj{t„) 



(18) 



where tn '= nh. 

The state space, M = ^^(Z^), is the space 
of Fourier coefficients; if a G M, then a = 
(^Ji' ^J2' ' ' where (ji, • • • ) is some ordering of 
Z^. Let (18) define the induced discrete-time evo- 
lution map, Th : M ^ M. At this point, we could 
exactly reduce this problem to example [2] by restrict- 
ing our attention to a finite-dimensional subspace of 
M. 

Let the observable space, J^, be the C- valued func- 
tions on M. A family of observables on M, param- 
eterized by X G B, is given by (16); namely, fixing 
we have for any a G M 



(19) 



Hence, the temperature at a point G B is a linear 
observable on the space of Fourier coefficients. 

Assume the temperature can be measured at a fi- 
nite number of points (iCi, iC2, • • • , ^k) in B and take 
the finite collection of observables {fx^ , • • • , /x^) 
fined by (19). Then the action of the Koopman op- 



erator on this set of observables is 



U 



where jUj = exp(- 



-ATT^c'WjWlh). 



(20) 



□ 



III. KOOPMAN MODE ANALYSIS 

A. Eigenfunctions and Koopman Modes of U 

Thus far, we have avoided putting structure on 
the function space J^. When is a vector space, 
the Koopman operator is linear. It, therefore, makes 
sense to study its spectral properties as this will give 
us insight into the dynamics of the system, simi- 
lar to the case of linear finite-dimensional systems. 
We make the further assumptions that is a Ba- 
nach space under some norm, ||-||, and that U is 
a bounded, and hence continuous, operator on this 
space. 

Let {01, . . . , 0n} be a set of eigenfunctions of [/, 
where n = 1,2,..., or 00, not necessarily forming a 
complete basis set for T. In the discrete-time case, 
we have that 



In the continuous-time case, the A's are eigenvalues 
of the generator U of the Koopman semigroup, {U^}. 
The eigencondition is then 



[U*^i]{p)=e^''Mp) 



(22) 



SO that {e'^*} are the eigenvalues for the Koopman 
semigroup. 

We first note two simple properties of eigenfunc- 
tions: their algebraic structure (Prop. |5| and their 
role in spectral equivalence of the systems (Prop. [t]). 

Proposition 5 (Algebraic structure of eigenfunc- 
tions under products). Assume T is a subset of all 
C-valued functions on M that forms a vector space 
which is closed under pointwise products of func- 
tions. Then, the set of eigenfunctions forms an 
Abelian semigroup under pointwise products of func- 
tions. In particular, , ^2 ^ eigenfunctions 
of U with eigenvalues Ai and X2, then 0i02 is an 
eigenfunction of U with eigenvalue A1A2. 

Furthermore, if p ^ and (j) is an eigenfunc- 
tion with eigenvalue X, then (jf is a eigenfunction 
with eigenvalue X^, where 0^(x) := (0(x))^. If (j) is 
an eigenfunction that vanishes nowhere and r G M; 
then (jf is an eigenfunction with eigenvalue X^ . The 
eigenfunctions that vanish nowhere form an Abelian 
group. 



Proof. Assume U(j)i = Ai^i and U(j)2 = X2(f 
put iIj{x) = 0i(x)02(^)- In discrete time. 



and 



[u^,]{p) = x,Up)' 



(21) 



[U^]{x) = ^(T{x)) = MT{x))MT{^)) 

= [U^i]{x) [UCI,2]{X) = XiX2M^)Mx) 

= AiA2V^(x). 

Hence, the set of eigenfunctions is closed under 
pointwise products. An analogous computation 
holds for continuous time. 

Note that constant functions are eigenfunctions 
at eigenvalue 1. Hence the constant function that is 
equal to 1 everywhere is an eigenfunction of U and 
acts as the identity element. Combining this with 
the above closure property and standard properties 
of pointwise products of functions shows that the set 
of eigenfunctions is an Abelian semigroup. 

Let = A(/) and fix p G M+. Then 

[ur]{x) = r{Tx) = {c^{Tx)Y = {x^{x)r 

= A^0^(x). 

If (f) vanishes nowhere, then (l)~^{x) := l/0(x) is well- 
defined. Then, the above chain of identities remains 
valid for r G M replacing p G M^. Hence, the Abelian 
semigroup of eigenfunctions also contains all of its 
inverses. Therefore, the set of eigenfunctions van- 
ishing nowhere is an Abelian group. □ 
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Example 6 (Analytic observables of sta- 
ble/unstable systems). Let x = Ax, with x,A G C 
and |A| ^ 1. Then ^^{x) = e^^x. Let = x. 

Then 

[/7V](^) = = (t){e^^x) 

which imphes that (j) is an eigenfunction of U . By 
proposition [sj any (j)n{x) := {(j){x))'^ = x'^ is an 
eigenfunction of U with eigenvalue A^. 

Let f{x) be an analytic function. Then f{x) = 

T^CnX"^ = X)cn</>n(^), whcrc = There- 
fore, 

[^/](^) = ^Cn[^</>n](:r) = ^ A-Cn(/>n(:r). 

□ 

The second property shows the spectral equiva- 
lence of topologically conjugate transformations. 

Proposition 7 (Spectral equivalence of topologi- 
cally conjugate systems). Let S : M ^ M and 
T : N ^ N be topologically conjugate; i.e., there 
exists a homeomorphism h : N ^ M such that 
S o h = h o T . If (j) is an eigenfunction of Us with 
eigenvalue \, then (j) o k is an eigenfunction Ut ctt 
eigenvalue A. 

Proof. Fix X G M and let y G N be such that x = 
h{y). The result follows from the chain of equalities: 

X{ct>oh){y)=X^{x) = [Us^]{x)=^[S{x)] 
= [UT{^oh)]{y). 

□ 

Example 8 (Topological conjugacy of diagonaliz- 
able systems). Let i/*^^^ = {y^\y^^y where the 
superscript {k) indexes time, and let = Ty^^\ 

where T is a matrix. Assume that T has eigenvectors 
Vi, V2 at eigenvalues Ai, A2 such that Vi 7^ e^, where 
Bj is the canonical basis vector. liV = [vi, ^2], then 



is an eigenfunction of [/a at eigenvalue 



after defining new coordinates x^^^ 
y-iy(fc)^ we get 



(fc) Jk).J 



,x 



Ai 0' 
A2 



x\ ' 



=: A 



(fe) 

Jb 9 



The maps A and T are topologically conjugate by 
AF-i = V-^T. 

Note, that (j)i{x^^^) = xf^ and (})2{x^^^) = 
are eigenfunctions of /7a at eigenvalues Ai and A2 
respectively. By proposition [5] and example [6) we 
have that 0m,n(ic^^^) := [(t)i{x^^^)Y'[(t)2{x^^^)T = 



A7^A2. By proposition j?! (l)m,n^V~^ is an eigenfunc- 
tion of Ut at eigenvalue A2^A2 , where V~'^ has taken 
the place of h in proposition j?! □ 

Now, assume / G J-* is an observable in the closed, 
linear span of a set of linearly independent eigen- 
functions (recall n could be finite or infinite). 
Then 



(23) 



for some constants Ci{f) G C. The dynamics of / 
are particularly simple: 

n 

[uf]{p) = mp)) = j2ci{f)MT{p)) 



= ^q(/)[C/</.,](p) 

i=l 
n 

= ^Kci{f)(l)i{p), 

and similarly 

n 



(24) 



(25) 



The extension to vector-valued observables F = 
(/15 • • • , /k)^7 where each fi is in the closed linear 
span of the eigenfunctions, is trivial: 



[u'f]{p) = Y,xTMp) 



C^{fl) 



Ci{fK)_ 



(26) 



where Q(F) := [q(/i), . . . q(/k)]^. Motivated by 
(26), we have the following definition. 



Definition 9. Let (f)i be an eigenfunction for the 
Koopman operator corresponding to the eigenvalue 
A^. Given a vector- valued observable F : M ^ 
the Koopman mode, Ci{F)^ corresponding to is 
the vector of the coefficients of the projection of F 
onto spanj^i}. 

Remark 10. The importance of defining Koopman 
modes with respect to eigenfunctions, rather than 
eigenvalues, becomes apparent when we consider 
vector-valued observables and non-simple eigenval- 
ues. For example, let A have a two-dimensional 
eigenspace and let and ^2 be a basis for it. 
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Let /i = ci(/)i +C2(/>2 and /2 = 03(^2 be scalar-valued. 
Define F = (/i,/2)^. The Koopman modes corre- 
sponding to and (j)2 are 



Ci{F) 



and 



C2{F) 



respectively. Note that both of these Koopman 
modes have A as the associated eigenvalue. There- 
fore, if the Koopman mode was defined with respect 
to the eigenvalue A, then it would not be a well- 
defined object. However, when the eigenspace is 
one-dimensional, there is no confusion in saying "the 
Koopman mode corresponding to A" . □ 

Remark 11. The definition of Koopman modes can 
be carried over with a slight modification to gen- 
eralized eigenfunctions. When an observable can 
be expanded in terms of only eigenfunctions, the 
Koopman modes are time-invariant objects. How- 
ever, when a generalized eigenfunction is present in 
the expansion, the Koopman modes become time- 
dependent objects. For example, let (j) be an eigen- 
function and ip a generalized eigenfunction of U cor- 
responding to A 7^ 0: 

/70 A0 and Uip = (j) ^ Xip. 

Let F = Ci{F)(j) + C2{F)\Ij be a vector- valued ob- 
servable. Then 

U^F = (Ci(F)A^ + C2{F)kX^-^) (j) + \^C2{F)iIj 
= X' (^Ci(F) + IC2{F)^ cP + \^C2{F)^ 

for k > 0. The Koopman mode for cj) at time k is 
the time-dependent quantity Ci{F) + kX~^C2{F). 

However, since Koopman modes of generalized 
eigenfunctions have not been treated in the litera- 
ture, whenever we refer to Koopman modes in this 
paper, we implicitly mean a Koopman mode corre- 
sponding to an eigenfunction. □ 

To be completely explicit, the eigenfunctions are 
C-valued observables on the state space M and the 
eigenpairs (Ai,^^) depend only upon the dynamics 
(M, T) and the function space J^, not on a partic- 
ular observable. The Ci(-)'s can be thought of as 
a mapping from the observable space into a vector 
space V; for example, in (26) above, Ci maps into 
C^. The map F 1— )► (j)iCi[F) is then a vector- valued 
projection operator onto the subspace spanj^i}. 

Remark 12. Given as a Banach space of scalar 
functions, one can ask for conditions on the geomet- 
ric multiplicity of A, i.e., dimension of the eigenspace 
^A- The general answer to this question depends on 
the dynamics T and on the particular space of ob- 
servables T chosen. We can give an answer for the 



case most studied in literature, when T preserves a 
measure /i with F = L'^{M^ii). In this case, all the 
eigenvalues of the associated Koopman operator U 
are on the unit circle. 

When T is an ergodic transformation (i.e., when 
any measurable set S invariant under T is either of 
zero or full measure), all eigenvalues of U are sim- 
ple (see Petersen, §2.^d3Sl). When T is not ergodic, 
the state space can be partitioned into ergodic sets: 
minimal invariant sets S such that the restriction 
T\s : S ^ S to any S is ergodic. Since all ergodic 
sets S are disjoint, they support mutually singular 
functions from As a result, the number of linearly 
independent eigenfunctions of U at any particular 
eigenvalue A is bounded from above by the number 
of ergodic sets in the state space. The number of 
such ergodic sets is highly dependent on the charac- 
ter of dynamics. The partition into ergodic sets, the 
ergodic partition, will be discussed in more detail in 
Section IVAl 

Note that the computational method discussed 
later in Section HIB 2| assumes fixing an initial con- 
dition pq G M, which effectively selects the ergodic 
set 5* C M which contains the point po- Then 
the Koopman operator L^l^ acting on Lp'{S^ ijl\s) has 
simple eigenvalues, where /ils- is the ergodic measure 
on the ergodic component S. 

There is an interesting relation between ergodic 
dynamics and Proposition |5] when the space of ob- 
servables = L'^{M^ fi) is defined with respect to an 
ergodic measure fi. Given an eigenfunction (j) e T 
with eigenvalue A, Proposition |5] guarantees that 
is an eigenfunction with eigenvalue A^, as long as 
(j)^ G T. If the eigenvalue is periodic (A^ = A 
for some /c > 2), then (j)^ lies in the eigenspace 
Ex. Ergodicity guarantees the simplicity of Ex^ 
and hence there exists a non-zero c G C, such that 
ll^^-c^ll^ =0, in L2(M,/i) norm. □ 

Example 13 (Linear systems^^). We first look at 
the case when the dynamics are given by a linear 
map, A : M ^ M, on some finite-dimensional, 
inner-product space M; i.e., Xm-\-i = Axm- Sup- 
pose A has a complete set of eigenvectors, de- 
noted by {vi, . . . , Vn}^ with corresponding eigenval- 
ues {Ai, . . . , A^}. Let {wj}i be the eigenvectors of 
the adjoint A* with eigenvalues {Xj}i^ normalized 
so that (vj^Wk) = Sjk- Consider the observable de- 



fined 



j{x) = (x^Wj). Then 



[U(l)j]{x) = (l)j{Ax) = (Ax.Wj) = {x,A*Wj) 

= i^^^j^j) = i^^^j) = >^j(l>j{x). 

(27) 

We see that is an eigenfunction of the Koopman 
operator. However, the functions do not ex- 

haust all of the eigenfunctions of U. For example. 
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by proposition |5j (l)j{x)(f)k{x) = (ic, Wj) (ic, Wk) is an 
eigenfunction of U. In particular, g{x) := (x^Wj)^ 
is an eigenfunction of U with eigenvalue for any 

keN. 

Let F be the vector- valued observable defined as 
F{x) = X when x G span{i;i, . . . , i;^}, where i < 
and zero otherwise; i.e., F acts as the identity on a 
subspace of M spanned by the first £ eigenvectors 
and has the complement of that subspace as its ker- 
nel. Then 

i £ 

F{x) = ^j) = E (^)^^- (28) 



and 



(29) 



From these expressions, we see that the eigenvec- 
tor, Vj^ of the linear map A is the Koopman mode, 
Cj{F)^ corresponding to □ 

Example 14 (Rotations of the circle). Let the state 
space be the interval M = [0, 1). Let uj G (0, 1) and 
define T : M ^ M by 



T{p)=p- 



mod 1. 



(30) 



Note that there is a natural identification of M with 
the circle T = R/27rZ and the functions on M with 
27r-periodic functions on R. It is well-known that 
if cj G Q, then every initial condition is periodic 
and if uj is irrational, then the trajectory starting 
from any initial condition densely fills M. Note that 
the dynamics preserve the Lebesgue measure. Let 
T = I/^(M), be the space of Lebesgue integrable C- 
valued functions on M and consider the observable 
Mp) = e^'^^'^'P, neZ. Then 



A27Tnuj 



Kip)- 



(31) 



Therefore, for any n ^ Ij^ (j)n is an eigenfunction of 
U with eigenvalue = e*^^^^. Since the trigono- 
metric polynomials are dense in L^(T), then for 
^(p) = E„6z^(")e''™PeLi(T), 

[Uh]{p) = /£(n)e*2-»-</,„(p). (32) 



If F = (/i, . . . , JkY is vector- valued observable, 
then 



[C/F](p) = ^e'2-"-,/,„(p) 



/i(n) 



\jK{n). 



(33) 



Hence the vectors of Fourier coefficients are the 
Koopman modes of the system. □ 



Example 15 (Partial differential equations). We 
continue with the heat equation example from above 
(example [4| . Recall that a map was defined on the 
space of Fourier coefficients by 



(r(a)),-=exp(-47rV||j||^/i)a,-. 



(34) 



Note that the canonical coordinate projections, 
^j(a) := ttj, are eigenfunctions for the Koopman 
operator, with eigenvalues A^ = exp(— 47r^c^ \\3\\a 
by the computation, 

TO(«) = 0i(r(a)) = (n«))i 

.2„2 



exp(-47r c ||j||2M«j 



(35) 



= exp(-47rV||j||^/i)0,-(a). 

For the observable fx{o) := X]jGZ2 we get 

[U'^U]{a)= ^ A7(/>,-(a)e^2-^'-. 

Suppose we can only measure the temperature at 
a finite number of locations {iCi, . . . ^Xk}^ then for 

F — {fxi 7 • • • 7 fxK ) 5 



[U^F]{a)=Yx^^j{a) 



oi27rj-£Ci 



(36) 



In the expression above, each Koopman mode, 

^i27rj-£Ci 



(37) 



is just a "shape" function on the physical space B = 
[-1/2,1/2] X [-1/2,1/2]. This stresses the point 
that the eigenfunctions are defined on the state space 
M while the Koopman modes are functions in the 
output space. □ 

The development thus far has only focused on the 
case when an observable is in the closed linear span 
of some set of eigenfunctions of the Koopman oper- 
ator. No assumption was made on whether this set 
was a complete set for U or even if U possessed a 
complete set of eigenfunctions. One could ask what 
conditions we could impose on the system that are 
sufficient to guarantee that U has a spectral decom- 
position. This is the case for measure-preserving dy- 
namical systems, as we now explain. 

Let ^ C M be the attractor of the dynamical sys- 
tem and fi the unique invariant measure supported 
on A. Often, /i will be a so-called physical measure. 
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These types of measures exhibit the important prop- 
erty 



1 r 



fe=0 



(38) 



for any continuous observable f : M ^ C and for 
Lebesgue-almost every p G M belonging to a pos- 
itive Lebesgue measure set ^ C M containing the 
attractor (see Youngl^ for a more detailed discus- 
sion). Such measures are important in applications 
since they guarantee the existence of well-defined 
time-averages even when an experiment starts with 
initial conditions not on the attractor. In such cases, 
we can restrict our attention to the dynamics and ob- 
servables on the attractor and recover all the asymp- 
totic behavior of the system. 

The situation on the attractor is quite nice when 
we consider the function space = {A^ ji) . The 
restriction of the dynamics to the the attractor, 
T\j^ : A ^ can be shown to be invertible /i-almost 
everywhere. The restriction of the Koopman opera- 
tor to the attractor, U\j\,: L'^{A^fi) /i), can 
then be defined by U IaI = / ^ ^U- this case, 
the operator is unitar}^^^^, implying that all of the 
eigenvalues lie on t he unit circle and the eigenfunc- 
tions are orthogona P^ * ^^ * ^^ 4 

Since U\j^ is unitary, there exists a spectral 
resolution^^ 



uUf 



SI 



XdE{\)f 



(39) 



where ^ is a project ion- valued Borel measure on the 
unit circle; i.e., for any Borel set 5* in the unit cir- 
cle, E{S) is a projection operator. The measure E 
is supported on the spectrum of /7|^. E can be de- 
composed into two measures, Ep and Ec^ that are 
supported on the point spectrum and the contin- 
uous part of the spectrum, respectively. For any 
/ G /i), the spectral resolution becomes 



(40) 



where Pj : L'^{A^ij) L'^{A^ijl) is the orthogonal 
projection onto the the eigenspace corresponding to 



the eigenvalue \j 



and Ec is the projection- 



valued measure corresponding to the continuous part 
of the spectrum. Either term on the right-hand 
side of (40) could be zero depending on whether 



the operator has no point spectrum or no contin- 
uous part of the spectrum. When the eigenvalues 
are simple, we get for a vector- valued observable 



Jo 



(41) 



where (j)j is the eigenfunction corresponding to the 



eigenvalue Xj 



and we have used PjF{jp) 



<pj{p)Cj{F) (when is non-simple, this identity 
does not necessarily hold). Finally, we note that the 
constant functions on the attractor are eigenfunc- 
tions ofU\A at eigenvalue 1 and that Pq/ •= Xa / 
defines a projection onto the constant functions. 
Then, (|4T| becomes 



[f/|$iF](p)= / F{p)d,,{p) 



J2TT0k 



dE,{0)F{p), 



(for further details on this decomposition consult 
Mezicl^. 

Therefore, for the case of measure-preserving 
transformations or those systems possessing a phys- 
ical measure, the asymptotic dynamics of the Koop- 
man operator are given by the contributions of three 
components: (1) the average value of the observable, 
(2) the portion admitting a Koopman mode expan- 
sion (the part corresponding to the point spectrum), 
and (3) the contribution of the continuous part of the 
spectrum. The third component, dEc{0)F{p)^ we re- 
fer to as the Koopman mode distribution, or the KM 
distribution for short. Whereas the Koopman mode 
expansion (the second component in the decompo- 
sition) is fairly well-understood with respect to its 
relation to the physics of a problem, the contribu- 
tion of the KM distribution does not enjoy the same 
level of understanding and has received little atten- 
tion, so far, in the literature. 



B. Computation of Koopman Modes 

In the examples that were presented thus far, it 
was fairly easy to determine the eigenpairs of the 
Koopman operator and the corresponding Koopman 
modes. For a general system, however, things will 
not be so easy. This section will discuss a few meth- 
ods to compute the projection of an observable onto 
the eigenspaces of the Koopman operator, from both 
the theoretical and numerical viewpoints. 
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Theoretical Results 



The first tool is given by the following theorem. 
It is a special case of a result found in Yosidaf^. 



Theorem 16. Let T be a Banach space and U : 
T ^ T . Assume \\U\\ < 1. Let A be an eigenvalue 
ofU such that |A| = 1. Let U = \~^U and define 

k=0 

Then Ak converges in the strong operator topol- 
ogy to the projection operator on the subspace of U- 
invariant function; i. e, onto the eigenspace Ex cor- 
responding to A. That is, for all f ^ 

1 ^"^ 

hm AKf = lim ^ E = ^a/. (43) 

/e=0 



Proof. Since each (pi is an eigenfunction, 
spanj^i, . . . , is a [/-invariant subspace, so 
the restriction of U to this subspace is a finite- 
dimensional linear operator and can be represented 
with a matrix. Any G spanj^i, . . . , 

can be written as fn = YlJLi(^j{fn)<l^j- Then 

fn - YliZl Ci{fn)(l)i e span{0^-,...,^^}. U 
restricted to spanj^j, . . . , (j)m} has eigenval- 
ues {Aj,...,A^}. Then Xj^U restricted to 



spanj^j, . . . , ( 



has eigenvalues {1, • • • , 



}• 



where Px : T ^ Ex is a projection operator. 
Proof. See Yosida^^ or Krengell^. 



The modulus of any element of this set is < 1. 
Taking the average, as in (43), using the operator 
XJ^U restricted to spanj^, . . . , gives the 

projection onto span of A~^ [/-invariant functions. 
The AJ^ [/-invariant functions are just elements of 
spanj^j}. Then 

k=0 \ i=l / 



□ 



Consider the case when the eigenvalues are simple 
and I All = ••• = |A^| = 1 and |An| < 1 for n > ^. 



is just Cj{fn)(t)j' This is equivalent to (45 ) when F = 
fn- The extension a vector- valued F is obvious. □ 



Then, A.- 



for some real cuj , when j < i. For 



vector- valued observables, the projections defined by 
(43) take the form 



Remark 18. Theorem 17 is a simple consequence 



K-l 



iCj{F) = lim 



-i27TUJi k 



[U^F], (44) 



of theorem 16 The case of F having elements in a 
generalized eigenspace is more difficult and is treated 
in forthcoming work by the authors. □ 

Remark 19. Analogous expressions hold for con- 
tinuous time, with (44) and (45) replaced by 



1 

for j = 1,...,^. Hence, Theorem [16] reduces to (j)jCj{F)^ lim — / e-''^''^^^[U^ F]dt. 

ht pvnppt for thosp T^oo / Jq 



Fourier analysis, as one might expect, for those 
eigenvalues on the unit circle and the projections 
can be computed with any implementation of a fast 
Fourier transform. For further discussion, consult 
Mezic and Banaszukl^ or Mezic^^. 

When an observable is a linear combination of a 
finite collection of eigenfunctions corresponding to 
simple eigenvalues, we get an extension of the above 
theorem to eigenvalues not having unit modulus. 

Theorem 17 (Generalized Laplace Analysis). Let 

{Ai, . . . , \m} be a (finite) set of simple eigenvalues 
for U ordered so that |Ai| > ••• > |A^| and let 
be an eigenfunction corresponding to A^. For 
each n G {1,...,A^}^ assume fn - M ^ C and 
fn G span{(/)i, . . . , (pni}- Define the vector-valued ob- 
servable F = (/i, . . .JnV- 

Then the Koopman modes for F can be computed 



(46) 



and 



= lim _ 



respectively. 



ri:' 



-Xjt 



i-1 

U'F-Y,e^''^^C^{F) 



dt. 

"(47) 
□ 



The first thing to note is that for theorem 17 and 
its continuous-time analogue a set of eigenvalues is 
needed; they are not computed as part of the theo- 
rem. To get the projections (pj Cj{F)^ the more un- 
stable modes must be subtracted off of the dynamics 
before the time-average is computed. 

The most common case for which Theorems [l6] 
(and their continuous-time analogues) are 



17 



via 



= lim ^ A. 

/c=0 



K-l 



J-1 

U^F-Y,^%C,{F) 



(45) 



and 

applied occurs when we restrict our attention to a 
compact invariant subset C of the basin of attraction 
for some attractor A (Z M and take for the prod- 
uct i^^(w4, ji) xl-L{C)^ where ja is the unique invariant 
measure supported on the attractor and H(C) is the 
space of analytic functions on C In this case, an 
eigenvalue for U satisfies |A| < 1. 
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Example 20 (Harmonic oscillator). Consider the 
Harmonic oscillator 

Pl{t)=P2{t) 
P2{t) = -Uj'^piit) 

Letting p{t) = {pi{t),p2{t))~^ , the solution flow is 



p{t) = $*(p(0)) 



COS out - sin out 
-uj sin ut cosut 



PiiO) 

P2{0) 



Note that this system is divergence- free and there- 
fore preserves volume in the state space; all eigen- 
values of U have modulus 1. 



Pi(0) 

P2{0) 

PiiO)- 



P2{0)/OJ 

-wpi(O) 



Let F{p{t)) = p{t). Equation (46) is nonzero only 
for LOi = uj and W2 = — w- Then for Ai = e*'^ and 
A2 = e~*", we get 

MpmciiF) 

and 

Mpm)C2{F) 



ou 



) 


"1" 




iou 



1 

2 

_ 1 

~ 2 

respectively. We can write 



Pi(0) 

P2{0) 

Pi{0)- 



+ i 



P2{0)/lo 
-wpi(O) 



.P2(0)^ 





1 


) 


—iou 



[U'F]{p{0))=F{p{t)) = <^\pm 



6i(p(0))Ci(F) 

"^V2(P(0))C2(F), 



or more explicitly, 



piit) 

P2{t) 



= e-^\(pm-^'-^ 



-iujt ' 



pi(0) + i 



P2(0) 



We recognize the familiar normal mode expansion 
for the harmonic oscillator. The normal modes 
are the Koopman modes Ci{F) = [1, ioj]^ and 
C2{F) = [1, —iou]^^ whereas the Koopman eigen- 
functions (/>i,2(p(0)) = §(Pi(0) TP2(0)/^) are the 
terms of the normal mode expansion that are func- 
tions of only the initial conditions. □ 

While in principle the projections onto the sta- 
ble and unstable modes can be computed directly 
using Theorem [iTj difliculties arise when we move 
past simple cases. If an explicit representation of U 
is not known for the observable F, we would need to 



compute the projections numerically. For |A| 7^ 1, 
this would require a numerical implementation of 
a Laplace transform; these are generally unstable 
computations, precluding the direct numerical im- 
plementation of Theorem [T7| when F has stable or 



unstable Koopman modes. Therefore, Theorem 17 
is more suited as an analytical tool, rather than a 
numerical one. 



2. A numerical algorithm: the Dynamic Mode 
Decomposition (DMD) 

Usually, we do not have access to an explicit rep- 
resentation of the Koopman operator. The behavior 
of the operator can only be ascertained by its action 
on an observable and usually at only a flnite num- 
ber of initial conditions. Thus, we are led to consider 
data-driven algorithms for computing the Koopman 
modes. By data, we mean a sequence of observa- 
tions of a vector- valued observable along a trajectory 
{T^p}. The following algorithms use these sequences 
of observations to approximate both the eigenvalues 
and the Koopman modes of U without having to 
numerically implement a Laplace transform. This is 
accomplished by flnding the best approximation of U 
on some flnite-dimensional subspace and computing 
eigenfunctions of the resulting flnite-dimensional lin- 
ear operator. The notion of best approximation will 
be made clear in what follows. 

Fix a vector-valued observable F : M ^ and 
consider the cyclic subspace 



/Coo = span{/7^F}g 



:0 5 



(48) 



that is, /Coo is the space of vector- valued observables 
in which flnite linear combinations of elements from 
{U^F}kLo are dense. 

Fix an r < 00 and consider the Krylov subspace 



JCr = span{/7^F}^ 



1 

k=0' 



(49) 



We will assume that {/7^F}^~q is a linearly inde- 
pendent set so that these functions form a basis 
for JCr. Note that UJCr C /Cr+i, so that, in gen- 
eral, JCr is not /7-invariant; it is only invariant if 
W^F e span{C/^F}^I^. 

Let Pr : J-^ JCr he di projection from the space 
of vector- valued observables onto JCr. Then 



PrU\)C^ : JCr ^ JCr 



(50) 



is a flnite-dimensional linear operator. This operator 
has a matrix representation, : ^ C^, in the 
{/7^F}^~Q-basis. Note that this matrix is dependent 
upon (1) the vector- valued observable, (2) the num- 
ber of time-steps r used (the dimension of the Krylov 
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subspace), and (3) the projection Pr used which is 
specified by the type of approximation we choose to 
make; in the following algorithms, this projection is 
the least-square approximation for evolution from a 
single point p G M. 

If (X^v) is an eigenpair for A^, where v = 
{vo, . . . , Vr-i^ G C^, then (j) = XI^Iq F] is an 
eigenfunction of PrU\)c^. By restricting our atten- 
tion to a fixed observable F and a Krylov subspace, 
we have reduced the problem of finding eigenvalues 
and Koopman modes of the Koopman operator to 
finding eigenvalues and eigenvectors for a matrix . 

A standard method for computing eigenvalues of 
a matrix is the Arnoldi algorithm and its variants. 
These are iterative methods relying on Krylov sub- 
spaces. The basic idea behind these algorithms is 
to project the matrix onto a lower-dimensional sub- 
space and to compute the eigenvalues of the lower- 
rank matrix. If the projection has a nice representa- 
tion, then the eigenvalue problem for this lower-rank 
approximation can be efficiently solved. 

The standard Arnoldi algorithm^ assumes we 
have a matrix A : whose eigenvalues and 

eigenvectors we want to compute. Starting from a 
random vector b G of unit norm, we form the 
Krylov subspace 

JCr := span{6, A6, . . . , A'^'^b}. 

Assuming full rank, an orthonormal basis {Qj}i for 
JCr can be found using a Gram-Schmidt procedure 
applied to {A-^6}J~q. Orthogonalization and renor- 
malization is usually performed at each step j. Let- 
ting be the matrix formed from the orthonormal 
basis, we get the relation 

H, = Q^AQ, 

where is of upper Hessenberg form and has the 
interpretation of the orthogonal projection of A onto 
JCr- can be diagonalized efficiently. The eigen- 
values of approximate the r eigenvalues of A of 
largest magnitude. Implementations of this algo- 
rithm use various additional methods to ensure nu- 
merical stability. 

By applying the Arnoldi algorithm, we implic- 
itly assume there exists a matrix A whose evolution 
A^6 G matches the evolution [U^F]{p) G for 
k — 0, . . . , r. Unfortunately, since we do not have 
an explicit representation of the Koopman operator, 
we cannot use the standard Arnoldi algorithm. This 
is due to the orthogonalization and renormalization 
performed at each step, which is equivalent to chang- 
ing the observable F at each step. If gr^. is the vector 
formed from normalizing the component of [/7^F] 
orthogonal to span{[/7-^F]}Q~^, then there is some 
G \ M ^ such that = G{p). Hence we are 



never looking at the action of the Koopman operator 
along a single trajectory and observable, precluding 
the use of the Arnoldi algorithm with data obtained 
from simulation. 

A variant of the Arnoldi algorithm, utilizing com- 
panion matrices, was ffrst described in Ruhe^^. The 
algorithm was popularized in the Fluids community 
by Rowley et al.^ and Schmid^. In Schmid^^, the 
algorithm was dubbed the Dynamic Mode Decom- 
position (DMD), whereas Rowley et al. related the 
algorithm to the approximation of Koopman modes. 

The strength of the DMD algorithm is that it only 
requires a sequence of vectors {6/e}^^Q, where 



bk := U'F{p) G e 



(51) 



for some ffxed F : M ^ and ffxed p e M. The 
algorithm gives the best approximation, at the point 
p G M, of the projection (pC{F) onto the eigenfunc- 
tion (eq. (26)). As will be seen, this corresponds 
to a speciffc choice of the projection operator Pr ap- 
pearing in ([50|. 

To derive the DMD algorithm, let 



[6o, . . . ,6r-l]- 



The columns of are the C^-vectors resulting from 
the point evaluations of the {/7^F}-basis for the 
Krylov subspace JCr at the point p G M. 

In general, br will not be in the span of the 
columns of K^. In this case, br = Yl^jZo ^j^j + Vr^ 
where the Cj's are chosen to minimize the C^-norm 
of the residual 77^. This corresponds to choosing 
the projection operator Pr appearing in (50) so that 



PrU^F is the least-squares approximation to U^F 
at the point p e M diS measured by the C^-norm; 
i.e.. 



\\[U^F]{p)-Pr[U^F]{p)\\^^ = 



r-l 



< 



r-l 



br -"^djbj 



for any other {do, • • • , <^r-i}- 

Since br = K^c + 77^, where c = (cq, . . . , c^-i)''', 
we get 

UKr = [61, ...,br] = [61, . . . , 6r-l, Kc + TJr], 

or equivalently 

UKr = KrK + Vr^^ (52) 
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where e = (0, . . . , 0, 1)"^ G and 



■0 
1 
1 







Co 
ci 

C2 



1 Cr 



(53) 



is the r X r companion matrix; it is the matrix rep- 
resentation of PrU in the {/7^F}^~Q-basis. 
Diagonahze the companion matrix: 

= V-^AV, (54) 

where A is the diagonal matrix of eigenvalues and 
the columns of are eigenvectors of A^. Inserting 
the expression for A^ into (52) and multiplying on 
the right by gives 



Define 



Then (55) becomes 



(55) 



(56) 



(57) 



For large enough m, it is hoped that ||?7^e'''V~-'^ || is 
small. If that is the case, then /7E ^ EA and the 
columns of E approximate some eigenvectors of U 
and the diagonal elements of A approximate some 
eigenvalues of U. Note that ||?7^e'''V~"'^ || = when- 
ever r > since then the columns of are linearly 
dependent which implies that 77^ = 0. 



be defined as in (54) 
column of 



Definition 21. Let A and E 

and ([56|, respectively. Let Wi be the z 
E and be the i^^ diagonal element of A. Then Wi 
is called an empirical Ritz vector and A^ is called an 
empirical Ritz value. 

Each empirical Ritz vector approximates 
4^i{p)Ci{F)^ the projection of F onto some eigenvec- 
tor and the empirical Ritz values approximate 
the corresponding eigenvalues of U. For this reason, 
we will generally refer the empirical Ritz vectors 
and values as the Koopman modes and eigenvalues 
computed by the DMD algorithm although this is 
not strictly true and loosens the terminology. 

Remark 22. The above algorithm is very much tied 
to the initial condition chosen. This dependence 
arises since the empirical Ritz values and vectors are 
formed using a Krylov subspace that is generated by 
a sequence of vector- valued observations along a fi- 
nite trajectory having initial condition p G M. A 



given initial condition may not reveal the full spec- 
trum and different initial conditions can reveal dif- 
ferent parts of the spectrum. For example consider 
the dynamical system. 



j Aip/c, Pk <0 



(58) 



where < Ai < 1 < A2. Let (j)i{p) = min{p, 0} and 
4^2{p) = max{0,p}. These are eigenfunctions of the 
Koopman operator at eigenvalues Ai and A2, respec- 
tively. Let F{p) = p = (t>i{p) + 02 (p)- Analytically, 
we can decompose the observable into a sum of pro- 
jections onto eigenspaces: F{p) = PiF{p) + P2^(p), 
where Pi is the projection onto span(/)^. Choosing 
an initial condition p < and applying the DMD al- 
gorithm only computes the projection onto the sta- 
ble mode; the DMD algorithm only reveals PiF{p). 
Similarly, choosing p > only reveals the unstable 
mode, P2F{p). Therefore, the DMD algorithm may 
only reveal a subset of the spectrum of the Koopman 
operator and the corresponding Koopman modes. 

It should also be remarked that if F ^ span{(/)^} 
for some eigenfunction 0^, then the DMD algorithm 
will not reveal that mode. This is often the case for 
natural choices for a set of observables, as was the 
case in the linear system example above (see example 
13)p.[9|. □ 

The version of the DMD algorithm described 
tends to be numerically ill-conditioned. The is due 
to A^6o converging to the eigenspaces correspond- 
ing to the largest magnitude eigenvalues, resulting in 
the columns of becoming nearly linearly depen- 
dent. A robust version of the algorithm has been de- 
scribed in Schmidts. It amounts to first computing 
a singular value decomposition (SVD) of and pro- 
jecting A onto the Krylov subspace using the SVD 
basis (for details, consult Schmid^^). Chen et alJ^ 
discusses variants of the Dynamic Mode Decomposi- 
tion, relates it to discrete Fourier transforms, and in- 
troduces an "optimized" DMD algorithm that com- 
putes an arbitrary number of modes from data. 



C. Applications of Koopman Modes 

The theory of Koopman modes has led to a num- 
ber applications in the literature. Broadly, the uses 
of Koopman modes can be classified under two head- 
ings: model reduction and coherency. Model reduc- 
tion deals with extracting the spatial features of just 
a few Koopman modes and attempting to under- 
stand the physics of the system just based on those, 
neglecting the details contained in other Koopman 
modes. The notion of coherency, on the other hand, 
deals with how observables relate dynamically with 
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respect to a Koopman mode. Coherency is always 
defined for a (not necessarily proper) subset of ob- 
servables and an eigenvalue A. The subset of observ- 
ables is coherent for A if the dynamics are identical. 
This reduces to checking if the initial magnitudes 
and phases of Cx{F) are the same for each observ- 
able in the subset. The following definition, first 
appearing in Susuki and Mezic^^, makes this pre- 
cise. 

Definition 23 (Coherency between Koopman 
Modes). Consider a vector- valued observable F : 
M C^, where F = (/i, . . . , fmV and fjiM^C 
for j = l,...,m. Let {Ci(F), . . . , Q(F)}, £ < oo, 
be a collection of Koopman modes of interest. Note 
that Ci{F) = {ci^u...,Ci^mV ^ C"^- Fix ei,e2 > 
and consider ji, j2 ^ {1, • • • , Then, fj^ and fj^ 
are (ei, €2)-coherent (with respect to the chosen 
Koopman modes) if 

(i) Ikijil - \cij^\ \ < ei, and 

(ii) \Zcij^ -^Cij^l < £2. 

for ah i = 1,...,^. □ 

The choice of two epsilon values in the above def- 
inition allows the practitioner to set the tolerances 
of the magnitude and phase independently. This is 
useful when one can tolerate more variation in either 
the modulus or the phase and still call the modes co- 
herent. Thus Cij^ is coherent with Q^j^, if the com- 
plex number Cij^ is contained inside some rectangle 



centered at Cj 



'',J2 



The examples we present are necessarily a subset 
of those that exist in the literature. 



1. Power system J^^J^ 

Koopman mode analysis has seen application in 
the analysis of power systems. In Susuki and 
Mezic^^^^, the authors used Koopman mode anal- 
ysis to identify coherency in the short-term swing 
dynamics of multi-machine power systems, with the 
New England Test System and IEEE Reliability 
Test System - 1996 being used as test cases for the 
methodolog}!^^!^. In Susuki and Mezicl^, the au- 
thors used Koopman mode analysis to identify pre- 
cursors to the so-called coherent swing instability 
of power systems where a group of generators syn- 
chronously loses coherency with the rest of the sys- 
tem after a local disturbance. We will focus only 
upon the identification of coherency, as pursued in 
Susuki and Mezic^^'^^, since it underlies the work 
in Susuki and Mezicf^ as well. 

The New England system is a 39-bus sys- 
tem having 10 synchronous generators, while the 



IEEE systems has 73 buses and 99 synchronous 
generator#^^'. As the analysis of the two systems 
is the same, we focus on the simpler New England 
system. 

The swing dynamics of the New England system 
were given by the following systems of differential 
equationi^^^ 



Hi duoj 
irfb dt 



Pmi — GiiEf 

y^^EjEjCjj, 



-DiUj, 

10 



(59) 



J = l 



where Cij is the coupling term 



Cij — Gij cos{Si 



5,) 



Bij sm{5i - Sj) 



In this model, i = 2, ... 10 indexed the generators, 
6i was the angular position of the rotor of genera- 
tor i relative to bus 1, and uji was the rotor speed 
of generator i relative to that of bus 1; Di was the 
damping coefficient of generator i, Ei was the volt- 
age of the generator, and Pmi the mechanical input 
power; Ga was the internal conductance of generator 



while Gi 



IBij was the transfer impedance 



between generators i and j; Hi was a pe r unit time 
inertia constant and a frequency^^^. The vari- 
ables Hi^ Ei, Di, fij and p ower loads were speci- 
fied during simulations^^lMl^ xhe variables Ga, Gij, 
and Bij were computed using power fiow computa- 



tions (see Susuki and Mezic^^^''^, and the references 
therein, for full numerical details). The state space 
for each generator {i = 2, . . . , 10) was the cylinder 
C = [— TT, tt] X R and the full state space M for this 
system was M = Cx---xC = C^. Each generator 
exhibited a stable equilibrium at ((^*,co'* = 0), for 
some (5*, computed using a power fiow computation. 

Let an observable fi : M ^ R be given by 
fi{6,uj) = uJi, where S = ((^2, • • • , ^lo) and uj 
(a;2, . . . ,^io)- The vector- valued observable chosen 
for Koopman mode analysis was F = (/2, . . . , /lo)^, 
so that 



F{S,u;) 



002 



^10 



This was a physically-relevant observable since in 
practice one meas ures the rotor speeds for each gen- 
eration planlPEl. 

The system was evolved for a short time period 
with a disturbance from the equilibrium state local- 
ized at rotor 8. Study of the resulting trajectories 
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showed that gene rators 2, 3, 6, and 7 were a co- 
herent groupP^l^ the four generators exhibited re- 
sponses in angular frequencies uji having the same 
amphtude and phase. 

The DMD algorithm was applied to the same sim- 
ulation data. The Koopman modes of interest were 
those that had the largest norms and correspond- 
ing growth rates |Aj|. Figure [l] shows a plot of 
the magnitude and phase of each component of the 
three Koopman modes with the largest growth rates 
and magnitudes. These were labeled as modes 7, 
8, and 9 and corresponded to frequencies of 1.3078, 
1.0962, and 0.3727 Hz, respectively. For each mode 
j = 7, 8, 9, there were amplitudes Aji and phases aji 
for the generators i = 2, . . . , 10. In the figure, the 
amplitudes and phases for modes 7, 8, 9 are plotted 
in plotted with symbols *, x, and o, respectively. 
Number labels within the plots specify the gener- 
ator. It is seen that generators 2, 3, 6, 7, and 9 
are coherent with respect to mode 8, while all but 
generator 9 are coherent with respect to mode 9. 
Therefore, generators 2, 3, 6, 7, and 9 were coherent 
with respect to both modes 8 and 9, as one found 
with visual inspection of the trajectories. While not 
done in this paper, the amplitudes and phases repre- 
sented in this way allow using a number of clustering 
algorithms to automatically identify coherency. 
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FIG. 1: Power Systems - amplitudes and phases for 

the components of the three Koopman modes 
having largest growth rates and norms. Koopman 
modes and eigenvalues were computed using the 
DMD algorithm. The largest modes are labeled as 
modes 7, 8, and 9. The amplitudes and phases of 
the components of mode 9 are plotted with the 
symbol o. The numbers in the plot correspond to 
component of the mode. Modes 7 and 8 are plotted 
similarly, expect with symbols * and x, 
respectively. Generators 2, 3, 6, 7, and 9 are 
coherent with respect to mode 8. All generators, 
except generator 8, are coherent with respect to 
mode 9. 



2. Jet in CrossflovPJ 

One of the earliest applications of Koopman 
modes was to the study of nonlinear fiuid fiows. 
Rowley et al.^^ introduced the concepts to the fiu- 
ids community and demonstrated the methodology 
by computing a subset of the Koopman modes, us- 
ing the DMD algorithm, for a jet in a crossfiow. 
The jet in a crossfiow configuration is a common 
way of mixing the jet fiuid with a uniform crossfiow. 
The crossfiow moves parallel to a fiat plate and the 
jet is injected through an orifice in the plate. It is 
known that such a system can exhibit self-sustained 
oscillation^^, and Koopman modes were used to au- 
tomatically identify the relevant frequencies and cor- 
responding three-dimensional fiow structures. 

The fiow field was studied in a (Z/^^, Z/^, Z/^) — 
(75,20,30)^0 computational box, where Sq was the 
displacement thickness at the crossfiow inlet. The 
incompressible Navier- Stokes equations over a fiat 
plate were solved using a Fourier- Chebyshev spectral 
method with a grid resolution of 256 x 201 x 144 (see 
the reference for all the simulation details). There- 
fore, each point in the state space M is a sequence 
of Fourier- Chebyshev coefiicients. The vector- valued 
observable, F : M ^ R^, was chosen as the veloc- 
ity measurements of the fiow field at the grid points. 



Hence, m = 3(256 x 201 x 144) ^ 2.2 x 10^ (three 
velocity components at each grid point). The situ- 
ation is similar to the heat equation example above 
(example |4j p. [6| where the state space was a space 
of Fourier coefiicients and the observable was a heat 
distribution on a square. 

Remark 24. The choice of state space M and a 
basis for the observables are not unique. For a fiuid 
experiment on some physical domain B, we assume 
that the solutions exists in some function space. 
Usually this is the space of finite energy fiows so that 
a velocity profile u{x) exists in L'^(B,dx). Depend- 
ing on the boundary conditions, many different bases 
for L^(B,(iic) exist. The particular basis chosen de- 
pends on the computational method used to solve 
the Navier- Stokes equations. Above, this method 
was a Fourier-Chebyshev spectral method; the re- 
sulting basis for L^(B, dx) consisted of trigonometric 
and Chebyshev polynomials. The state space M was 
the sequence space of Fourier-Chebyshev coefficients 
corresponding to L'^(B^dx) functions. However, if 
periodic boundary conditions for the ffuid ffow are 
assumed, the trigonometric polynomials could be 
used to represent L'^(B^dx) and M would be the 
space of Fourier coefficients. Therefore, the choice 
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of solution method for the Navier- Stokes equations 
determines the basis functions and implicitly defines 
the state space M that represents the system. □ 

An initial velocity profile and boundary condi- 
tions were specified for simulations and this con- 
figuration corresponded to a fixed p e M. Leting 
bk := [U^F]{p), the DMD algorithm was performed 
on the sequence of observables {6200, ^202, • • • , ^700 }• 
The reason for the delay in the start time of the se- 
quence was to neglect the transient terms. 

The top row of figure [2] shows the time signals of 
the streamwise velocity recorded by a sensor probe 
close to the wall, just downstream from the jet ori- 
fice and a probe located downstream and on the jet 
shear layer. In the bottom row, the spectral content 
of the time signals are shown in black, whereas the 
spectral peaks identified by the DMD algorithm are 
shown in red. Note that for the probe near the wall 
(left column of figure [2| the signal contains only low- 
frequency components, whereas the signal near the 
jet contains both low- and high-frequencies. The 
probes are local measurements and only pick up a 
subset of the full spectrum for the fiuid fiow. The 
Koopman modes are global objects and, as shown 
in the figure, the DMD algorithm identifies both 
the low- and high- frequency components of the fiow 
field. 

Figure |3] shows the spectrum of the Koopman op- 
erator as computed by the DMD algorithm. Most 
eigenvalues lie on the unit circle implying that the 
fiow field is near an attractor. The time-averaged 
flow (steady-state component) corresponds to A = 1 
and is indicated in blue in the left image of figure 
[Sj The rest of the eigenvalues have colors smoothly 
varying from red to white with the colors corre- 
sponding to the magnitude of the associated global 
mode. Red corresponds to large magnitudes for the 
Koopman modes, white to low magnitude. The mag- 
nitudes are given by the total energy of the mode (2- 
norm). The right image of figure [s] shows the mag- 
nitudes of the Koopman modes at each frequency. 
The color scheme is the same as for the left image. 

If we order the Koopman modes in order of de- 
creasing magnitude, mode 1 corresponds to the 
time- average fiow and the rest come in complex- 
conjugate pairs; modes 2 and 3 correspond to 
complex- conjugate eigenvalues and have the same 
magnitude. Figure [4j shows the streamwise velocity 
components of mode 2 (left) and mode 6 (right). 
Each mode oscillates at a single frequency, with 
mode 2 corresponding to a high-frequency {St = 
0.141) and mode 6 to a low- frequency {St = 0.0175); 
St is the Strouhal number. These correspond to the 
tallest red line and the left most red line in the bot- 
tom row of figure [2j respectively. In both images 
of figure |4j the red surfaces correspond to positive 



streamwise velocities and blue surfaces to negative 
streamwise velocities. Mode 2 is associated with 
shear layer vortices with additional vortices extend- 
ing toward the wall. Mode 6 has large structures 
along the wall associated with shedding of the wall 
vortices. This shedding of wall vortices is coupled to 
the main jet body as indicated by the mode having 
structure along the jet body. 



3. Self-sustained oscillations in a turbulent cavit^^ 



Seena and Sungl^ investigated the causes of self- 
sustained pressure oscillations in fiuid fiow over a 
cavity (see figure |6] for a profile of the experiment 
geometry) by using the Dynamic Mode Decompo- 
sition algorithm. We denote this domain by B. 
The goals of the study were to identify the vorti- 
cal structures that drove the hydrodynamic oscilla- 
tions and obtain dynamical information about those 
structureP^. Both thin and thick incoming bound- 
ary layers were studied (Reynold's numbers at the 
cavity of 12000 and 3000, respectively); only the tur- 
bulent case {Re = 12000) is covered here. 

As with the case of the jet in a crossfiow example, 
the fiuid evolution was governed by the incompress- 
ible Navier- Stokes equations corresponding to the 
specified cavity geometry (see Seena and Sung for 
details on the computational domain). Abstractly, 
the state space M was a sequence space of coeffi- 
cients for basis functions on B. Since the system 
was not periodic in all dimensions, the basis func- 
tions were not the 3-dimensional trigonometric poly- 
nomials (see Remark 24). A family of observables, 
parameterized by points in B, was given by functions 
mapping a point in M (a sequence of coefficients) to 
the fiuid pressure at a point x eM. 

The solutions of the Navier- Stokes equations were 
computed using the Crank-Nicolson method using 
a second-order central difference scheme in space^. 
Cavity ffows at the high Reynold's number were 
simulated with a Large Eddy Simulation (LES)^. 
Broadly speaking, LES filters the governing equa- 
tions by removing the small-scale structures and 
replacing them with models. Consult Seena and 
Sung^^ for full details of the simulation parameters 
and Germano et al.^ for a discussion of LES. 

The vector-valued observable chosen for analysis 
was the fiuid pressure at each of the computational 
grid points in and above the cavity. The DMD al- 
gorithm was applied to a sequence of 124 fiow snap- 
shots recorded after allowing transients to decay^. 
Figure [5] shows the Koopman eigenvalues resulting 
from the DMD algorithm. On the left, almost all 
of the eigenvalues are seen to be on the unit circle, 
implying that the fiow field is on or near an attrac- 
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FIG. 2: Jet in Crossflow - The top row is the time signal of the streamwise velocity for a probe near the 
wall, downstream of the jet orifice (left) and for a probe downstream, in the jet trajectory (right). In the 
bottom row, the spectral content of the corresponding probe signals are shown in black. In red, the part of 
the spectrum of the Koopman operator captured by DMD algorithm. Only the positive frequencies are 
shown since eigenvalues occur in complex-conjugate pairs. (Original in Rowley et al. Journal of fluid 
mechanics by Cambridge University Press. Reproduced with permission of Cambridge University Press in 
the format reprint in a journal via Copyright Clearance Center.) 
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FIG. 3: Jet in Crossflow - (left) The spectrum of Koopman operator as identified with the DMD algorithm. 

Most of the Koopman eigenvalues are on the unit circle. The mode corresponding to the time-averaged 

flow (corresponding to A = 1) is indicated in blue. The other eigenvalues are colored from red to white 
based on the total energy of the associated Koopman mode. Red corresponds to high-energy modes, white 

to low-energy modes, (right) The magnitudes of the Koopman modes at the each frequency. The color 
scheme is the same as for the image on the left. (Original in Rowley et al.S^, Journal of fluid mechanics by 

Cambridge University Press. Reproduced with permission of Cambridge University Press in the format 

reprint in a journal via Copyright Clearance Center.) 



tor. Colors in the plot correspond to the total en- 
ergy of the corresponding Koopman mode. On the 
right, the total energies of the Koopman modes at 
each frequency are shown. The four dominant peaks 
are labeled. Only positive frequencies are marked 
since eigenvalues occur in complex-conjugate pairs. 
The mode labeled as 1 occurs at cj = and corre- 
sponds to the steady component of the flow. The 



corresponding Koopman mode is shown in figure 
|6ja). Solid, black lines correspond to high-pressure 
regions, whereas broken, gray lines correspond to 
low-pressure regions. Koopman modes 2, 3, and 4, 
corresponding to the labeled peaks in figure jsjb), 
are shown in figure |6jb)-(d). Modes 2 and 3 corre- 
spond to a; = 4.6 rad/s and 3.5 rad/s^ respectively. 
The exact frequency corresponding to mode 4 was 
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FIG. 4: Jet in Crossflow - Koopman modes 2 (left) and 6 (right) corresponding to high- {St2 = 0.141) and 
low- {Ste = 0.0175) frequencies, respectively. St is the Strouhal number. Red contours correspond to 
positive streamwise velocities and blue contours to negative streamwise velocities. Reproduced from 

(Original in Rowley et al.^^, Journal of fluid mechanics by Cambridge University Press. Reproduced with 
permission of Cambridge University Press in the format reprint in a journal via Copyright Clearance 

Center.) 



not reported, but was less than 10 rad/s. The low- 
and high-pressure regions of the modes oscillate at 
a fixed frequency and suggest self-sustained oscilla- 
tions in the cavit}!^. These results were consistent 
with results reported in literature^^. The study suc- 
cessively used the DMD algorithm to identify the 
large vortices in the cavity driving the self-sustained 
oscillations. 

Remark 25. Figures |3] and |5] show that the Koop- 
man spectrum, as computed by the DMD algorithm, 
has some eigenvalues inside the unit circle, even 
though an initial block of data points was discarded 
before computation so that transients could die out. 
They are most likely spurious eigenvalues resulting 
from the finite truncation of the data and the DMD 
algorithm and not due to slowly decaying modes con- 
tained in the data. Unfortunately, no detailed in- 
vestigation has been done on such eigenvalues; most 
attention has been focused on the modes correspond- 
ing to eigenvalues on the unit circle and those inside 
the unit circle have been ignored. □ 



4. Energy efficiency in buildingJ^^^ 

Energy use and efficiency in buildings has re- 
ceived much attention in recent years and Koopman 
mode analysis has recently seen application in this 
field, as well. A researcher is generally interested 
in temperature distributions in buildings, as this 
can tell much about HVAC and controller perfor- 
mance or if the building is operating near its design 



pointP^E^. Such meas urements can also be used for 
model validatiorP^E^. 

Analysis requires understanding of the heat fiow 
in the building subject to forcing from weather or 
even human traffic between areas of the building. In 
the most general case, heat transfer equations must 
be solved for a complicated domain with varying 
boundary conditions. The interesting observables 
are the temperatures at each point in the building, 
B. The temperature distributions in the building are 
assumed to exist in some function space (L^(B,(iic), 
for example). A temperature distribution can be 
represented in a basis 62, • • • } for the function 
space. The state space M is the sequence space of co- 
efficients for these basis functions; i.e., a point p e M 
is p = (ci, C2, . . . ), for some coefficients G C. An 
observable maps a point in M to the temperature at 
a point X eM. When collecting real data, the tem- 
perature is measured at a finite number of points, 
usually dictated by the placement of the tempera- 
ture sensors by the building designers. The vector- 
valued observable of interest is the vector of the tem- 
peratures recorded by the building sensors. 

In the most general case, the full model for the 
building is too complicated to solve. A variety of 
simplifying assumptions and a numerical package is 
needed to compute solutions to the building system. 
EnergyPluJII' is a simulation package for modeling 
energy and water use in buildings. It is a free pro- 
gram offered by the U.S. Department of Energy and 
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FIG. 5: Turbulent cavity flow - (left) Koopman eigenvalues computed using the DMD algorithm for 
turbulent flow inside the cavity at Re = 12000. The eigenvalues are colored based on the energy of the 
corresponding Koopman mode, (right) The energy of the Koopman mode at each frequency, uj {rad/s). 
(Original in Seena and Sung^^, International journal of heat and fluid flow by Institution of Mechanical 
Engineers (Great Britain) Reproduced with permission of Institution of Mechanical Engineers] in the 
format reuse in a journal/magazine via Copyright Clearance Center.) 
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FIG. 6: DMD modes inside the cavity at Re = 12000 for modes labeled 1,2,3, and 4 in flgurejs] (Original 
in Seena and Sung^^, International journal of heat and fluid flow by Institution of Mechanical Engineers 
(Great Britain) Reproduced with permission of Institution of Mechanical Engineers] in the format reuse in 

a journal/magazine via Copyright Clearance Center.) 



is widely usedQ At a high level, an EnergyPlus 
building model consists of specifying the location of 
every surface in the building. A list of rooms is then 
created and the interaction of the rooms (between 
shared surfaces) is specifled. For example, heat con- 
duction through a certain type of material may be 
specifled for one wall, whereas radiative heat transfer 
is important for windows. Models can also include 
weather data, HVAC usage, room occupancy, and 



http : //appsl . eere . energy . gov/buildings/energyplus/ 



building water usage. The software makes the as- 
sumption that rooms are well-mixed so that no tem- 
perature gradients exist in the room. Under the as- 
sumptions of the software, heat transfer in the build- 
ing is approximated by a system of coupled ordinary 
differential equations. 

In Eisenhower et al.l^, Koopman mode analy- 
sis was used to decompose the temperature evolu- 
tion in a building into purely periodic global modes. 
The building investigated was the Y2E2 building 
at Stanford University and was modeled and sim- 
ulated in the EnergyPlus software package. The 
physical building has 2370 HVAC sensors recording 
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data throughout the buildin^^. The vector-valued 
observable chosen was the vector consisting of tem- 
perature readings taken at the sensors on the second 
fiooi^. Koopman mode analysis was used to vali- 
date the EnergyPlus model by comparing the most 
energetic Koopman modes for the simulation and 
the raw data. 

Figure [7| shows the spectral content of the Koop- 
man modes computed using the DMD algorithm 
with the sensor data on top and the EnergyPlus data 
on the bottom. The sensors are along the vertical 
axis with the period of the Koopman mode along 
the horizontal axis. A vertical streak in figure [7| cor- 
responds to a "global" mode; i.e., most sensors are 
affected by the Koopman mode at that frequency. 
In both images, a strong streak is seen at the 24 
hour period which is due to the daily forcing from 
the weather. Horizontal blue lines indicate near con- 
stant temperatures at those sensors and correspond 
to rooms served by their own fan units^^. 

Figure |8] shows the 24 hour Koopman mode for 
both the sensor data and the EnergyPlus model. 
The top row corresponds to the sensor data, the 
bottom row to the simulation. The magnitude of 
the mode is shown on the left with the phase on the 
right. Note that the magnitude and phase of the 
Koopman mode is reported in reference to the out- 
side air temperature (OAT). The relative magnitude 
of a mode is given in decibels (dB) 



|i^M,| =201og 
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Coat{F) 
and the phase is given in degrees 

AKM, = AC,{F)-ACoat{F), 

where Ci{F) and Coat{F) are the Koopman modes 
as computed by the DMD algorithm. For the sensor 
data, the magnitude of the 24 hour mode relative to 
the external temperature was about —6 dB through- 
out the building. This implied that the peak tem- 
perature oscillation inside the building was smaller 
than the magnitude of the temperature fiuctuations 
outside. For the EnergyPlus mode, the temperature 
magnitude inside was much closer to the external 
fiuctuations since the relative magnitude of the mode 
was about dB. Significant deviations in the rela- 
tive phases were also seen between the sensor data 
and the model. These discrepancies implied that the 
parameters used in the simulation model were incor- 
rect. The analysis lead to suggestions on the model 
parameters to modif}'^^. 

In practice, even under the simplifying assump- 
tions utilized in programs such as EnergyPlus, a de- 
tailed model of the building can be prohibitively ex- 
pensive to simulate over the time scales needed (e.g.. 



hours up to years). Methods for model reduction are 
in order. Zoning approximations are one approach. 
In a detailed EnergyPlus model, each room is treated 
as a unique thermal zone. Each room is assumed 
well-mixed so that the thermal properties are uni- 
form in the room. The idea of zoning is to lump ad- 
jacent rooms together in such a way that the thermal 
properties can then be assumed uniform across those 
rooms. This procedure results in less regions that 
need to be simulated in the model. Usually, zoning 
approximations are performed heuristicall}'^^. How- 
ever, model accuracy is quite sensitive to the zoning 
approximation usecP^. 



Georgescu et al.l^ used the notion of coherency 
between Koopman modes (def. 23 above) to create 
zoning approximations for buildings and studied the 
Engineering Sciences Building (ESB) at the Univer- 
sity of California, Santa Barbara as a test case for 
the methodology. The particular observable chosen 
was F = (/i, . . . , fmY 1 where each fjiM^ R, 
j = 1, . . . , m, represented the temperature of a room 
in the ESB. In the detailed EnergyPlus model, there 
were m = 191 zones. An additional physical as- 
sumption was imposed so that rooms were grouped 
into a zone if they were both on the same fioor and 
adjacent as well as e-coherent. 



The modes of interest were those having periods of 
one year, 24 hours, 12 hours, 8 hours, and 6 hours. 
These corresponded to the most energetic modes. 
Figure [9] shows the Koopman modes of the ESB 
EnergyPlus simulation for the three most energetic 
modes. The magnitudes of the Koopman modes are 
given in the top row. Units are degrees Celsius rel- 
ative to the average temperature in the room. The 
average temperatures were not reported. The bot- 
tom row corresponds to the phase in radians of the 
Koopman mode. Consider the block of four rooms 
in the center of the right hand side of the building 
(the ones that are dark blue in the magnitude plot 
of the year long mode). By visual inspection, the 
middle two rooms on this block would be lumped 
into a single zone. While the colors of those two 
rooms vary between each of the six images in fig- 
ure [9] within the same image those two rooms have 
the same color, and hence almost identical values. 
These two rooms satisfy the e-coherency definition 
for some small e > 0. 
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FIG. 7: Y2E2 building decomposition - Koopman 
spectra for the sensor data (top) and the 
EnergyPlus model (bottom). A large spectral 
content is seen at the 24 hour period due to 
outdoor forcing conditions. Horizontal blue lines 
correspond to near constant temperatures at those 
sensors. These sensors correspond to rooms that 
are served by their own fan units allowing a tight 
control of temperature. Reproduced from 
Eisenhower et alJ^. 




FIG. 8: Y2E2 building decomposition - Koopman 
mode for the 24 hour period, (top row) The 
magnitude and phase, respectively, of the raw 
data's Koopman mode, (bottom row) The 
magnitude and phase of the EnergyPlus model's 
Koopman mode. The magnitude units are given in 
decibels and phase in degrees. Both are relative to 

the external 24 hour mode. The discrepancy 
between scales shows the model mismatch with the 
real data. Reproduced from Eisenhower et al.-^. 
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FIG. 9: Zoning approximations of buildings - Koopman mode magnitude and phase for the modes having 
the largest magnitudes. Only one floor of the ESB is shown, (top) the magnitudes of the Koopman modes. 
Units are in degrees Celsius relative to the mean temperature of the room, (bottom) the phase of the 
Koopman mode in radians. Reproduced from Georgescu et al.l^ 



IV. ANALYSIS OF STATE SPACES USING 
EIGENQUOTIENTS 

Previous sections discussed the Koopman eigen- 
functions in the context of spectral decomposition 
of the Koopman operator, and we did not pay much 
attention to values that eigenfunctions take on the 
state space. In this section, we demonstrate that 
eigenfunctions carry information about transport 
between parts of the state space and provide a way 
to identify invariant sets. Going further, we will en- 
dow the collections of invariant sets with their own 
metric topology, which allows for associating smaller 
invariant sets into larger coherent structures. 

The material in t his section ap peared originally 
in a string of paperJ^EISHMIlTlig ^-^j^ techni- 
cal det ails g iven also in dissertations of two of the 
authors. 1^1^ 



A. State-space analysis of Koopman eigenfunctions 

The level sets of Koopman eigenfunctions at eigen- 
values |A| = 1 form invariant sets, and periodic and 
wandering chains of sets in the state space. The 
eigenfunction level-set partitions depend not only on 
the choice of the eigenspace Ex from which eigen- 
functions are taken, but also on the choice of a 
particular function in the eigenspace. However, an 



exhaustive procedure for analysis of level-set parti- 
tions can be devised, based on partition products, 
such that the result in the limit does not depend 
on the particular choice of eigenfunctions studied, 
but rather only on the eigenvalue A, which deter- 
mines the dynamics of level sets, e.g., periodic or 
invariant. In this section, we first focus on the er- 
godic partition, which corresponds to analysis of the 
invariant eigenspace of the Koopman operator, i.e., 
A = 1. The generalization to "periodic" eigenspaces, 
i.e., those corresponding to A = e*^^^ for cj G Q, is 
straightforward, and we present it at the end of this 
section. 

From this point onward, we will assume the mea- 
surable setup of the Koopman operator: T : M ^ 
M will be a measurable map between measure spaces 
(M, 05), with at least one invariant measure /i. The 
observables are, at the very least, a subset of 
complex- valued measurable functions. Given an 
invariant function ^, i.e., a function in invariant 
eigenspace of the Koopman operator (j) e Ei{U)^ we 
can form its level set partition ({<p) := {Sz G 
C}, where the level sets are Sz := {x G M : (j){x) = 
z}. Since is a measurable function, its level sets Sz 
are measurable sets. Additionally, they are invariant 
sets, by virtue that each Sz collects all state-space 
points on which <p takes the value z, i.e., if it con- 
tains X, it will also contain T(x), T^(x), . . ., since 
(p is constant along trajectories. A partition of the 
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state space into invariant sets is called a stationary 
partition. 

Unless the system is ergodic, the choice of <p G Ei 
is not unique (see Remark 12): often there will be at 
least two functions 0, G , which are not linearly 
related and yield two different partitions C(^), ClV^)- 
The product of two partitions CW^CH^)^ which con- 
tains intersections of sets in each of the partitions 
(see Fig. [lO|, is again a stationary partition. The 
product is finer than either of the original partitions, 
as sets in ({(j)) or C('0) can be recovered as unions of 
sets in ({(j)) V C('0)- Naturally, given a finite set of 
functions ^ = 1? • • • ? ^5 the product V/c=i Ci^'k) 
of the associated partitions can be formed incremen- 
tally. 



C(0) 



C(0)vCW 



M 









FIG. 10: Partitions of the state space M into 
({ip) and their product ({(/)) V CW- 

Given a partition ( of the space and some set A, 
the quotient map tt^ : M ^ A is any map that sepa- 
rates the partition (: 7r^(x) = 7r^(^) iff 35 G C such 
that x^y G S. In other words, all the points in any 
5* G C will be mapped to a single point p G A. The 
image set ^ := 7r(M) is termed the quotient ^ of the 
partition (. If we are given just the map tt^, we could 
reconstruct the partition ( from level sets of 7r<^, in 
terms of values of the quotient map codomain A. 

For a product partition = VfeLi C(</^/c)7 the 
most straightforward way to construct a quotient 
map is by arranging eigenfunctions in a vector- 
valued function to define ttk - M C^, TTxix) := 
(01 (x), ^2(^)7 • • • , 0k(^))- Notice, however, that the 
quotient map is not unique: any bijection h : 
A, e.g., a coordinate transform, can be used to form 
another quotient map through composition h o ttk- 
While the construction of ttk might appear aca- 
demic, it is significant in that it translates the prob- 
lem of constructing partition from the topologi- 
cal setting, where intersections of sets were used to 
check whether two points x, y are in the same prod- 
uct set, to an analytical setting, where the equality 
7^k{x) = TTxiy) can be checked, either numerically 
or analytically. We will make use of this fact heavily 
in Section ITbI 

Evaluation of eigenfunctions 0/^, which feature in 
construction of quotient maps, is easy when we are 
given a method for evaluating the projection of the 
space of observables onto an eigenspace Px : T ^ 
Ex. In this case, we can compute a dense set of 



eigenfunctions in Ex by evaluating Pxf for a dense 
set fk G T. Moreover, when is a Hilbert space, we 
can compute the projection on an orthogonal basis 
for J^, obtaining a (not necessarily orthogonal) basis 
for Ex. A practical way to numerically evaluate Px 
for |A| = 1 using trajectory averages is explained in 
Section IVCl 

For A = 1, projecting functions fk using Pi re- 
sults in an infinity of invariant functions that we 
can attempt to use in forming incremental level-set 
partitions. It is not trivial that the partition inter- 
section can be extended to the case of an infinity of 
U -invariant functions 6^. The first result of this sort 
was obtained by Sine^H for continuous observables 
T = C{M). Such setting is fairly restrictive for de- 
terministic dynamics. Therefore, we focus on a more 
general resuliP^ that holds for T = L^{M,ii). 
Assume we are given a countable, bounded family 
G £^1 C L^{M^ /i), for k — 1,2,..., whose span is 
dense in the invariant eigenspace E\ of the Koopman 
operator U . The incremental partition into level sets 



exists as the level set partition of the map TTg : M ^ 
^ defined as TTg := (^1,^2, • • • where is the 
space of bounded sequences. The map TTg is termed 
the ergodic quotient map^ with ergodic quotient •= 
7Te{C) as its image set. 

While the concept of the ergo dic partition was 
known since at least Rokhlin,^^ this computation- 
ally feasible construction is fairly recent. The 
justification for the name of the ergodic partition 
comes from analyzing restrictions of T to invariant 
sets S e Ce- Each S carries a measure /i^ such 
that the restricted dynamics T : S ^ S is ergodic 
with respect to /i^. By one of the several equivalent 
definitions of ergodicity,^^ this means that for any 

feL\S,^is) 



L 



/MS= lim (60) 

at almost every x G 5, with respect to jas- The 
measures fis are called ergodic measures and are 
extreme points of the space of invariant measures 
(by the Ergodic Decomposition Theorem). More- 
over, integrals against ergodic measures define the 
projection onto the invariant eigenspace Pif{x) = 
Is{x) f^l^s{x)^ where S{x) denotes the ergodic set 
containing point x. This is a consequence of the 
equality (60) along with Mean Ergodic Theorems,!^ 
e.g., von Neumann theorem for = L'^{M,/j.) or 
Yosida ergodic theorem for general Banach spaces, 
used in earlier sections as Theorem [TBI 
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The ergodic partition is unique up to /i-measure 
zero sets, i.e., given two stationary partitions ("i and 
(2 that satisfy the above properties, for any G Ci 
there exists S2 ^ (2^ such that jj.{SiAS2) = 0. Con- 
sequently, the partition ("e does not depend on the 
choice of functions (j)k used to construct it: differ- 
ent choices of functions just yield different represen- 
tations of the quotient map TTg, however, the limit 
partition will be the same. 

When there is only one T- invariant measure /i, up 
to multiplication by a scalar, the system is uniquely 
ergodic. In that case, the eigenspace Ei is one- 
dimensional, containing only a.e. constant functions. 
In turn, the ergodic partition contains a single set, 
and the ergodic quotient map maps almost all points 
(with respect to /i), into a single point in the se- 
quence space i"^. While uniquely ergodic systems 
might be of interest in general, they are not the tar- 
get of this approach. The ergodic quotient analysis 
is intended to help the analysis of state spaces which 
contain a lot of ergodic measures, e.g., systems with 
families of periodic or quasiperiodic orbits and sys- 
tems possessing interspersed regions of regular and 
irregular dynamics. 

For most conservative regular systems, ergodic 
partitions are similar to partitions of the state space 
into orbits. However, in zones of irregular dynam- 
ics, e.g., on strange attractors and in zones of strong 
mixing, it is difficult to gain intuition about be- 
havior of the systems simply from orbits. Orbits 
are inherently zero-dimensional (for maps) or one- 
dimensional sets (for flows), but in chaotic regimes, 
every one of them densely fills a bigger set, possibly 
even a set of a positive invariant measure, e.g., posi- 
tive volume. Moreover, in a mixing zone any two tra- 
jectories look nothing alike, yet they are contained in 
the same set and have the same statistical behavior. 
Through such reasoning, we might be interested not 
in description of trajectories themselves, but rather 
in description of minimal invariant sets containing 
individual trajectories. Precisely, the ergodic parti- 
tion C is a measurable hull of the decomposition of 
the space into orbits T^(x), i.e., a partition of the 
space into minimal measurable invariant sets that 
contain orbits. In this sense, ( is the appropriate 
counterpart of the state space portrait in the context 
of measurable dynamical systems. 

The understanding of the ergodic partition, de- 
rived from functions in the invariant eigenspace Ei^ 
can be extended to eigenspaces for which |A| = 1. 
When the eigenfunctions are taken from a "peri- 
odic" eigenspace for an eigenvalue with prop- 
erty = 1, for some 6 G N (period), the extension 
is straightforward. Level sets of (j) G Ex are then pe- 
riodic sets, which are arranged in b- periodic chains, 
{5„T(5,),T2(5,),...,T^-1(5J}. Similarly, if A = 



^i27TUj such that Q, the chains of level sets 
extend into infinity, forming wandering chains. The 
partition ({(j)) is again an invariant partition, in the 
sense that for each S G C{4^)^ ^('S') ^ CW- However, 
({(f)) is not stationary, as S are not invariant sets; it is 
an invariant 6-periodic partition or, when a; Q, an 
invariant wandering partition. The partition prod- 
ucts also generalize: if G Ex, ({(j)) V ("(V^) will 
also be a 6-periodic/wandering invariant partition. 

The invariant measures can be generalized to the 
concept of complex eigenmeasures,^^ which satisfy 

^^{T-'S) = A-V(5). 

Integrating against extrema of eigenmeasures en- 
ables us to evaluate Pa away from A = 1, and con- 
struct eigenquotient maps 7Tx{x) = (. . . , Pxfk-, - - ')i 
analogously to constructions of ergodic quotient 
maps from a basis fk for J^. In this sense, the ergodic 
quotient map is the eigenquotient map at A = 1. 

The eigenquotient maps collect basis functions for 
the eigenspace Ex- Note that this basis set might be 
overdetermined, since we do not know ahead of time 
what dim Ex is. In Remark \V2\ we have indicated 
that the dimension of Ex is bounded by the num- 
ber of mutually singular components of measure fi 
which was used to define the space of observables as 
T = L^(M, /i). These components are precisely the 
ergodic measures, and there are as many of them as 
there are ergodic sets. 

Given an eigenfunction (j) G ^^a, for |A| = 
1, its modulus and complex angle functions 
are respectively |^|, Z0, for which ^(x) = 
|(/)(x)| exp[z27rZ0(x)]. The modulus function is 
an invariant function, since U = oT\ = \X(j)\ = 
101 . On the other hand, e Ex is a. factor map: 
it conjugates the dynamics of T with a dynamical 
system -\-uj, evolving on a circle S^ , w ith uj as 
the angle of the eigenvalue A = re*^^'^,'^^El] since 

Z0(Tx) = Z0(x) +0;. 

To illustrate what these functions might reveal 
about dynamics, we will take the Chirikov Standard 
Map 



Xn+i = Xn^Pn + esin(27rXn) 
Pn+1 =Pn + esin(27rxn) 



(61) 



with e = 0.15. The system preserves the area mea- 
sure dm , therefore, the Koopman operator U : 
L'^{M,dm) L'^{M,dm) is unitary and its eigen- 
values all lie on the unit circle. Figure 11 shows 
level sets of an invariant eigenfunction, and modulus 
and angle of an eigenfunction at e^'^^^ for uj = 1/3, 
demonstrating the distinction between (j) ^ Ei and 
the modulus of G £^a- 
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(a) Invariant eigenfunction 

(j) e El. 
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(b) Modulus IV^I of a 3-periodic 
eigenfunction -0. Modulus of any 
eigenfunction for |A| = 1 is an 
invariant function. 
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FIG. 11: Level sets of eigenfunctions (j) ^ Ei and il) ^ Ex for A 



(c) Angle of a 3-periodic 
eigenfunction ijj. Angle of an 
eigenfunction at is a factor map 
for dynamics. (Only regions where 
angle is well defined, i.e., '0 ^ 0, 
are colored.) 

at = 1/3 for the Chirikov 



Standard Map at e = 0.15. (Functions were evaluated by computing Fourier averages, see Section IV C) 



Even though both |?/^| and (j) are invariant func- 
tions, IV^I is non-zero only over trajectories that are 
periodic with period 3. At e = 0.15, such trajecto- 
ries are mostly concentrated in two large groups of 
periodic chains near the middle of the state space. 
On the other hand, invariant functions that are 
not formed by taking absolute values of periodic 
functions do not necessarily make a distinction be- 
tween periodic and non-periodic dynamics. The an- 
gle function /^ip can be used to infer the order in 
which sets are ordered in the associated periodic 
chains. Even such simple visualizations can be use- 
ful aids in quick assessment of a new dynamical sys- 
tem, by identifying regions that are not dynamically 
connected. 



B. Geometry of eigenquotients 

In previous section, we explained how eigenfunc- 
tions of the Koopman operator connect to invariant 
and periodic structures in the state spaces of dy- 
namical systems. We have also indicated that the 
finest invariant partitions ^ can be represented in 
sequence spaces using quotient maps tt : M ^ £^ 
and their image sets, eigenquotients ^ = 7r(M). In 
this section, we will be treating geometry of sets ^, 
explaining it for the case of the ergodic quotient, 
which corresponds to invariant sets. As explained 
in the previous section, a generalization to periodic 
sets is straightforward. 

Identification of ergodic sets can be thought of as 
checking whether for two state space points pi , P2 ^ 
M it holds that 7r(pi) = 7r(p2) or not; essentially, we 
are using a discrete topology on ergodic quotient ^ 



to compare points in the state space M. We can, 
however, use other metric topologies on the ergodic 
quotient to extract additional information about the 
state space, e.g., to obtain a low-dimensional repre- 
sentation of dynamics, or identify functions acting as 
integrals of motion over invariant sets, not necessar- 
ily entire state spaces. By a low-dimensional repre- 
sentation of dynamics, we mean that we are looking 
for a minimal set of directions in the state space, in 
which we can move a point across the boundary of 
an ergodic set S and land in another ergodic set S' 
where trajectories look similar, on average, to those 
in S. 

As a motivation for the analysis of ergodic quo- 
tient we take the Morse theory analysis of Hamilto- 
nian systems. Consider two state spaces sketched in 



the top row of Figure [12] Each state space contains 
a region where trajectories are tightly layered to- 
gether, and where they can be parametrized using a 
single continuous parameter: distance from elliptical 
fixed points. For Hamiltonian systems, these prop- 
erties can be formally stated by inducing a topol- 
ogy on the state space based on the energy function. 
In other words, the "similarity" of neighboring er- 
godic sets that we are after is here represented by 
the similar values that the energy function attains 
on neighboring ergodic sets. 

Level sets of any monotonic function of the Hamil- 
tonian can be represented by Reeh graphs sketched 
next to the state spaces. Reeb graphs are a topolog- 
ical tool, featured in Morse theory, which analyzes 
manifolds through level sets of differentiable func- 
tions on them. The Reeb graphs provide a concise 
description of changes in topology of such level sets 
as the level, i.e., value of the function, is changed. 
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FIG. 12: Sketches of two Hamiltonian oscillators: a 
harmonic oscillator and a double-well oscillator. 
Top row: state space portrait. Bottom row: 
graphic representation of Fomenko-Reeb graphs for 
level-sets of Hamiltonian functions. Dashed arrows 
indicate how moving an initial condition across the 
level sets reflects on the Reeb graph. (Based on 
original in Budisic and Mezic^, Physica D: 
Nonlinear phenomena by North-Holland. 
Reproduced with permission of North-Holland in 
the format reuse in a journal/magazine via 
Copyright Clearance Center.) 



In dynamical systems, they have been used to study 
level sets of integrals of motion for integrable Hamil- 
tonian systems;^ in this context, they are known 
as Fomenko-Reeb graphs. Since integrals of mo- 
tion are invariants of projections Pi onto the A = 1 
eigenspace, in certain settings we could draw paral- 
lels between Fomenko-Reeb graphs and the ergodic 
quotient ^. 

The vertices in Fomenko-Reeb graphs correspond 
to those values of the Hamiltonian : M ^ R at 
which the topology of level sets of H changes. The 
edges connecting the vertices correspond to families 
of connected components of level sets that are ho- 
motopically related with respect to continuous vari- 
ation of values of H. Figure [12] illustrates what 
the Fomenko-Reeb graphs look like for two simple 
Hamiltonian systems. Without going into details 
of their construction here, an important feature of 
these graphs is that continuous segments in them 
correspond to one-parameter families of periodic or- 
bits. For the double-well oscillator, strands merge at 
the value of energy where two separate wells merge 
across the separatrix to form the outer family of pe- 
riodic orbits. Therefore, by looking at number of 
independent parameters in a Fomenko-Reeb graph, 
we can deduce the number of families of periodic 
orbits in the state space. 

Currently, Fomenko-Reeb graph analysis is formu- 



lated only for integrable Hamiltonian systems .^^ By 
analyzing the ergodic quotient, we can establish an 
approach similar in spirit, even for systems that do 
not have an explicit energy function, but might con- 
tain families of invariant sets in the state space. The 
ergodic quotient map collects, in a sense, all possible 
invariant functions, that can be thought as general- 
izations of energy. Unlike the case of Hamiltonian 
systems, where topology of the energy codomain was 
used to define topology on the state space, the er- 
godic quotient ^ is a subset of an infinite-dimensional 
sequence space which can be endowed with 
many euclidean-like topologies that are not all mu- 
tually equivalent. For this reason, we need to choose 
a metric structure (^, (i), where (i is a distance func- 
tion, so that we can obtain a context in which a gen- 
eralization of Reeb analysis of Hamiltonian systems 
to non-Hamiltonian systems is possible. 

The role played by the Hamiltonian function in 
Fomenko-Reeb graphs is taken up by the quotient 
map TT. The image space of Hamiltonian was R with 
its natural metric topology; the topology on the 
ergodic quotient ^ is too sensitive to establish anal- 
ogous arguments. Instead, we will use a Sobolev 
space metric which requires interpretation of the im- 
age under the quotient map 7r(x) as a spatial Fourier 
transform of ergodic measures jix- Within this set- 
ting, we seek to extract coherent onion layers in 
the state space by identifying connected components 
in The challenge lies in extracting a partic- 

ularly appropriate low-dimensional parametrization 
of ^ even though it is embedded within the infinite- 
dimensional space 

First, as observables which will be projected onto 
El we choose the normalized harmonic basis for J^, 
e.g., on M ^ [0,1]^ 

where k G and k ■ x = Yl'j=i ^j^j- This results 
in the representation of the ergodic quotient map tt 

^(p) = (..., [PiA](p),...) 
= (...,Ap(^),--0 

that maps the points of the state space p to Fourier 
coefficients fip{k) := Jj^ fkdjj^p of the associated er- 
godic measure fip. This opens up opportunities to 
use metrics available on the space of Fourier coef- 
ficients, e.g., weighted metrics, to compare dy- 
namics in different ergodic sets, and generalize the 
Fomenko-Reeb analysis from Hamiltonian systems 
to a more general class of systems. Note, however, 
that when we are interested only in identifying er- 
godic partitions using discrete topology on ^, i.e., 
ask whether x^y G M are elements of the same er- 
godic set 5* G C by checking if 7t{x) = 7T{y) , any 
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continuous basis of observables will be sufficient.^ 
The choice of basis is largely a matter of convenience 
of constructing other metric structures on ^, alter- 
native to the discrete metric. 

Due to boundedness of Fourier coefficients, the er- 
godic measures are elements of the Sobolev space 
whose norm can be defined as a weighted 
euclidean norm 



E 



iA(fc)r 



[1 + (2/!c7r)2 



(62) 



with the index determined as s = {D -\- l)/2 where 
D = dimM. We use as the symbol for the as- 
sociated Fourier coefficient space with the weighted 
euclidean norm. An excellent introduction to the 
Sobolev space theory is the classic textbook by 
Adams and Fournier^, which contains all the ma- 
terial relevant for this analysis.^ 

To establish a metric structure on the ergodic quo- 



tient, we use the 



\2-s 



-induced metric, and analyze 



continuity of ergodic quotient maps tt in it. The 
^ ]^2,-s structure induces a distance- 

like function on the state space 



(63) 



This function is not a true distance but a pseudo- 
distance precisely because the points at zero distance 
from each other are those that lie in the same ergodic 
set. Clearly, when the system is uniquely ergodic 
with respect to /i, the function dg evaluates to zero 
for /i-almost any pair of points pi,p2- 

Remark 26. The choice of the Sobolev space W^'~^ 
has a justification in comparison of ergodic measures 
interpreted as measures of residence times of trajec- 
to ries in measurable sets of a compact state space 
M.I^^Let B{x, r) denote euclidean balls in M and 
Xx,r their characteristic (indicator) functions. The 
quantity iip[B{x^r)] is the residence time of trajec- 
tory T^{p) in B{x,r), due to ergodicity of the sys- 
tem T : S ^ i.e., equality of integrals and time 



averages (60): 



1 ^"^ 

[B{x,r)] = / Xx,rdiiip = lim — V Xx^roT'^ip) 

Q N^oo iV ^ — ' 



where S C M is the element of the ergodic partition 
containing the initial condition p. 

The distance between trajectories originating at 
Pi •, P2 can then be formulated by integrating the dif- 
ference in residence times over all spherical sets, with 
R chosen such that M C B{x^R): 

d{pi,P2)'^ '= I / \lip^[B{x,r)]- iip^[B{x,r)\\^ dxdr. 
Jo Jm 

(64) 



The choice for the index of Sobolev norm s = 
{D -\- l)/2 results in the equivalence of norms, i.e., 
existence of constants a, /3 > such that for any 
Pi,P2 ^ M 

d{pi,P2) <ds = ll/ipi ^ ^ d{pi,p2). 



This argument shows that (64) is again only a 
pseudo-distance, like dg. We will return to this inter- 



pretation in Section |V| the formulation ( 62 ) is more 
useful in analysis of the ergodic quotient ^. 

The space (^, ||.||2 _g) can be analyzed computa- 
tionally, allowing for extraction of coherent struc- 
tures in state spaces of dynamical systems. We are 
looking to extract connected components of ^, i.e., 
curves C C In our motivational Hamiltonian 

examples, regions mapping to curves in Fomenko- 
Reeb graphs (see Fig. 12 ) contained an onion layer- 
ing of trajectories, which could be interpreted as re- 
gions of uniform dynamical behavior. We can visual- 
ize such regions by coloring all the points x G 7r~^(C) 
using the same color. 

In a recent paper, ^ we have presented an algorithm 
that performed the analysis described above. Using 
averaging of harmonic observables along trajectories 
(see Sec. IV C ), it evaluates an approximation of the 
ergodic quotient map tt on a set of initial conditions 
{Pn}n=i covering the region in state space that is 
to be analyzed. The connected components are ex- 
tracted based on pairwise evaluations of the induced 
Sobolev metric d{pi,p2) ^ -/^P2ll2,-s- 

The Diffusion Maps algorithnPISl was used to pro- 
vide a coordinate change for the ergodic quotient ^. 
The particular ergodic quotient map tt used earlier 
was constructed using averaged harmonic functions 
as coordinates. This choice was driven by desire to 
efficiently evaluate the Sobolev metric as a weighted 
euclidean metric. However, the harmonic basis set 
is chosen without any regard for intrinsic geometry 
of ^, and therefore it might not be the most efficient 
way of analyzing geometry of ^ . 

The Diffusion Maps algorithm interprets ^ as a 
heat-conductive object, and computes the modes of 
heat spread along it, resulting in a coordinate change 
^ : ^ ^ The f' distance between points ^(^) 
corresponds to the diffusion distance: an intrinsic, 
coordinate-independent distance along the ergodic 
quotient. This refiects the fact that the Diffusion 
Maps algorithm obtains intrinsic geometry of ^, re- 
gardless of the coordinate system that ^ was origi- 
nally represented in. Components of ^ are functions 
V^/e : ^ ^ R, ordered with /c G N according to the 
spatial scales over which they vary. For example, if 
^ is a simple line segment, as it is the case for the 
harmonic oscillator, il)k will just be Legendre poly- 
nomials: solutions of the heat equation with no-fiux 
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boundary condition. For more general ^, the diffu- 
sion modes ipk become more complicated, however, 
they retain the scale-ordering of harmonic functions. 
The ipk are good candidates for the analysis similar 
to the Reeb graph analysis for manifolds, which was 
previously noticed in medical image analysis.!^ 

As mentioned earlier, a bijection composed with 
the ergodic quotient map tt results in another ergodic 
quotient map. Therefore, by forming x •= ^ o tt, 
with components Xfc = V^fc ^ we obtain the er- 
godic quotient represented in a "good" coordinate 
set, which can be efficiently truncated to obtain 
low-dimensional representations of the ergodic quo- 
tient. The euclidean distance over a truncated set 
of diffusion coordinates approaches the diffusion dis- 
tance limit exponentially fast, in number of diffusion 
modes retainecP^. Since euclidean distance is the 
most common distance used in applied problems, a 
host of off-the-shelf algorithms can be used to post- 
process the ergodic quotient, e.g., a /c-means clus- 
tering or proximity graph analysis for extraction of 
connected components. 

We use the double- well potential to illustrate this 
analysis, with results presented in Figure [l3bj The 
Hamiltonian function for this system, which serves 
as the basis for the well-known Duffing oscillator, is 
given by 



1 



1. 



(65) 



in this analysis we chose k = 1^ b = 2. From Figure 



13b it is evident that the computations using met- 
ric II.II2 _s retain the desired intuition established by 



Figure [12] the diffusion coordinate X2 approximates 
the energy function of the system, while coordinate 
Xi discriminates between wells of the potential. The 
gaps in the numerical result are due to finite number 
of initial conditions x used to evaluate the ergodic 
quotient map, and irregularities are due to the finite 
averaging process used to evaluate the projection 
Pi; the Diffusion Maps algorithm can be adaptively 
tuned to tolerate such errors. When the state space 
is visualized using pseudo-coloring based on diffu- 
sion coordinates, the regions of uniform dynamical 
behavior are clearly distinguished. 

The presented process goes beyond Hamiltonian 
flows: as an illustration we use a periodically forced 
3d fluid flow based on the classical Hill's vortex 
flow.^ It is a solution to an ODE system on x = 
{R,z,0) C M+ X R X T: 
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where parameters c and e are swirl and perturba- 
tion strengths, respectively. Figure |13c| shows two 
vortices that exist at c = e = 0.3495 colored based 
on the colors assigned to components of ergodic quo- 



tient shown in Figure 13d 



These two examples demonstrate that functions 
Xk ^ ^ are invariant eigenfunctions for the Koop- 
man operator, but whose level sets resemble energy 
functions. In this sense, instead of finding a func- 
tion that would provide a "good" labeling of ergodic 
sets, as we did with energy functions in motivational 
Hamiltonian examples, we constructed a set of such 
a function from data generated by analyzing a basis 
for the invariant eigenspace of the Koopman oper- 
ator. The entire process is made computationally 
feasible by truncating the set Xki therefore discard- 
ing the fine-scale features in favor of comput ability, 
while retaining coarse-grained features in an orga- 
nized manner. 



C. Infinite-time averages as projections onto 
eigenspaces 

In our presentation of analyses of the state space 
using eigenfunctions of the Koopman operator, we 
have assumed we can compute eigenfunctions (j)^^^ 
by projecting an observable / onto the eigenspace 
at A using the projection operator Px. When A = 
g*27ru;^ i.e., for eigenvalues on the unit circle, we can 
evaluate the associated projection operator Px using 
infinite-time averages. For this reason, we restrict 
ourselves to the eigenspaces Ex for which |A| = 1 in 
this section. 

The focus of our interest will be the averages 



N ^ ' 

n=0 



7oT"(a;), 



(67) 



which we would want to extend pointwise into the 
limit as ^ 00. To describe the set of observables 
T which have well defined limits, we start by assum- 
ing that the map T conserves a measure \i on the 
Borel algebra ^ in M. For /i-integrable observables 
/ the finite averages of the form (67) can be ex- 



tended into TV 00 to result in well defined Fourier 
averages j^^^'. 

Theorem 27 (Wiener, Wintner). Let T : M ^ 

M preserve a measure fi on a measurable space M . 
Then the set of points T^{f) C M on which the limit 



~f^-\x) 



1 ^"^ 
lim — e" 

n=0 



VoT»(x) (68) 



(66) 



is well defined can be chosen for all f G V-^M^ji) 
and G independently of uo, and such that 
fi{M/^) = 0. 
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(a) State space of a double-well potential. Color is (b) Embedding of ^ for a double-well potential into 
the first diffusion coordinate xi? corresponding to first two diffusion coordinates. Colors are va lues of 

[b] xi for easier comparison witli[a| (Cf. Fig. 12 

center) 
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(c) Two vortices in the Poincare section of the 
periodic 3D Hill's vortex flow extracted based on 
continuous segments in ^. 



(d) Embedding of ^ for the periodic 3D Hill's vortex 
flow into first two diffusion coordinates. 



FIG. 13: Geometric analysis of the ergodic quotient for a planar Hamiltonian system with Hamiltonian 
(65) and a periodically forced 3d fluid-like flow defined by (66). (Original panels ([c| and ([d| in Budisic and 
Mezic^, Physica D: Nonlinear phenomena by North-Holland. Reproduced with permission of 
North-Holland in the format reuse in a journal/magazine via Copyright Clearance Center.) 



Proof. The original proof by Wiener and Wintner^ 
was found to contain an error. Several correct proofs 
have been compiled by Assani^^. □ 

Moreover, when L^(M, /i) contains a dense count- 
able set, the convergence set H can be chosen inde- 
pendently of fW^ When frequency a; = is chosen, 
Fourier averages (68) are referred to as ergodic av- 
erages: 



N-l 

f(x) := lim — V /oT^(x). 

n=0 



(69) 



As mentioned before, the concept of eigenmea- 
sures fi'p^ for map T provides means for evaluat- 
ing the projection operator, to which we now add a 



trajectory- wise formulation: 

JM 

A recent monograph by Wichtrey^^ provides a de- 
tailed analysis of existence of Fourier averages and 
their applications in linear, nonlinear, and control 
systems. 

While this presentation would perhaps suffice in 
the abstract, in practical settings we should take 
additional interest in: (a) the influence of choice 
observable / on the information contained in eigen- 
functions Pa/, (b) the "size" of the convergence set 
E on which the averages converge, and (c) the rate 
or error in approximating f^^^ using finite averages. 
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The influence of uj on the value of f^^\x) is easily 
seen by flxing / and x and studying the sequence of 
(complex) numbers f ^ T^(x). The average 



N-l 
n=0 



is then just a Discrete Fourier Transform (DFT) of 
the sequence a^. 

From the basic knowledge of Fourier analysis, the 
values /^^^ will be non-zero only for those uj that 
correspond to temporal modes present in the se- 
quence a^, which represent the trace of the observ- 
able / o T^{x). Whether contains a mode at uj 
depends not only on trajectory T^{x), but also on 
the choice of the observable. This is easiest to see 
when / is chosen as a constant function (if constant 
functions are in J^). In that case, 

fi^) = for all 

non-zero cj, and f'^^^ = /, regardless of the underly- 
ing dynamics T. 

Remark 28. Several concepts in this paper are as- 
sociated with Joseph Fourier's name. This is under- 
standable, as his ideas are at the root of all of them, 
but unfortunate as it might imply direct connections 
where there are none. 

Several recent work^^^^^ use the name harmonic 
averages for f^^\ in reference to harmonic analysis 
of dynamical flows. To avoid confusion with har- 
monic means of numbers, we decided to use Fourier 
averages instead, as do Wiener, Wintner, and As- 
sani. This name is connected with the temporal 
Fourier transforms of sequences, as explained above, 
since for flxed x e M and / G J^, the function 
UJ f^^\x) is the Discrete Fourier Transform of 
the sequence f o T^{x). 

Only through ergodicity, i.e., equality of space 
and time averages Jj^ f dfi = /, do we obtain a 
connection with spatial Fourier transform. Choos- 
ing fk{x) as harmonic functions e*^^^"^, where k is 
now a wave number, we can interpret any sequence 
k fk{x) in the ergodic quotient ^ as the spatial 
Fourier transform of the averaging ergodic measure 
jiix- We refer to the image set of the spatial Fourier 
transform as the spatial Fourier coefficients. 

The eigenquotient maps, i.e., x ^ 
(. . . , Pxf{x), . . . ), could be justiflably named 
Fourier quotient maps, in analogy to ergodic/Fourier 
average dichotomy. Despite our earlier choices of 
terminology,^ we expect this connection to hold only 
for eigenvalues A = e*^^*^, when Pxfk{x) = f^'^\x) 
As one could conceivably formulate an eigenquotient 
map for |A| ^ 1, in this paper we chose to use the 
term eigenquotient instead. 

Finally, the harmonic functions //c(x) = e*^^^'^, 
k e i.e., solutions to Laplace equation A/ = 



on a torus T^, are sometimes called Fourier func- 
tions or Fourier harmonics. We use these functions 
as observables, i.e., functions that are acted on by 
the Koopman operator, to facilitate evaluation of 
Sobolev norms. Out of all, this connection in name 
is the least signiflcant of all presented here, never- 
theless, we mention it in this remark to clarify our 
terminology. 

Remark 29. Starting from the same observable 
/, two different eigenfunctions at = can be 
computed as (f){x) := /^^^(x) and for any a; 7^ 0, 



ure 
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An example can be found in Fig- 
where the flrst two images were computed by 
Fourier averages of /(x,p) = sin(7rx — 7r/4) cos(67rp) 
at = and at a; = 1/3, after which the modu- 
lus was taken. Both observables are invariant with 
respect to dynamics. 

For a generic observable /, the main difference 



is that 



(x) is sensitive to cj-periodic dynam- 
ics: even though the information about the phase of 
cj-periodic sets is removed by taking the modulus, 

f^^"^ {x) will be zero on all points that are not in 

any of cj-resonant chains, which can be easily ex- 
plained by DFT interpretation. There is no such 
restriction on /^^\ whose value can vary between 
different invariant sets (see Fig. 11a). Conversely, 

2ero everywhere except on the period-3 is- 
land, where trajectories contain uj = 1/3 frequency 
components (Fig. |llb] ). 

The UJ values corresponding to the e*^^^ that are 
eigenvalues of the Koopman operator will result in 
non-zero eigenfunction /^'^^ for at least some observ- 
able / . The UJ values that are not frequencies of any 
of the eigenvalues will result in /^^^ =0. The fre- 
quency = 0, however, corresponds to eigenvalue 
A = 1 of [/, which is always in the spectrum of the 
Koopman operator. This frequency is of a particular 
interest, as the eigenfunctions /^^^ = / are invariant 
functions for the system. 

To illustrate the difference in eigenfunctions ob- 
tained by starting from different observables, we 
again used the Chirikov Standard Map (61) at e = 
0.15 . The averages were computed for = 0, which 
projects observables to the invariant eigenspace of 



U. As Figure [14] shows, the detail that level sets of 
averaged observables reveal about the state space is 
highly dependent on the starting observable, how- 
ever, when a "good" observable is chosen, it can 
reveal a lot about the state space. The lack of in- 
tuition about how to select such a "good" starting 
observable led to development of the ergodic quo- 



tient analysis, presented in Section [IV B| The diffu- 
sion coordinates Xk used therein can be interpreted 
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(b) Level sets of 
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FIG. 14: Three different observables fj (top row) and associated invariant functions fj (bottom row), 
obtained by ergodic averages along trajectories of the Chirikov Standard Map (61) for e = 0.15. The 

pseudocolor is the value of each function. 



as constructed observables that reveal detailed in- 
formation about the state space. 

The convergence set S was established to be of full 



measure ji by Theorem 27 The measure /i is a mea- 
sure preserved by the system, which is also used to 
define the space of integrable observables L^(M, /i) 
for which the Wiener- Wintner theorem holds. De- 
pending on the analyzed system, the existence of 
an invariant measure can be inferred through differ- 
ent arguments. For example, a system x = F{x) 
governed by divergence-free vector field V • F = 
conserves the volume- measure on the state space. 
In dissipative dynamics, the systems may conserve 
Sinai- Ruelle-Bowen measures on chaotic sets.^^ Fi- 
nally, on a compact state space M the only necessary 
assumption for existence of an invariant ji is conti- 
nuity of T, by the theorem of Krylov-Bogolyubov.l^ 

The Wiener- Wintner theorem holds for any in- 
variant measure ji. Therefore, by choosing the 
measure used to formulate the space of observables 
F = L^{M^fi)^ we influence the amount of informa- 
tion we can collect about the system by evaluating 
Fourier averages and the size of the set S. In ap- 
plied contexts, we would often want to include as 
many open sets as possible in the support of mea- 



sure fi chosen: for volume-preserving systems, vol- 
ume of the state space is often a good choice. Cau- 
tion is still needed, as there could be a /i-zero, yet 
dense, non-convergence set S^. Nevertheless, in our 
experience, simulations of dynamical systems that 
model physical phenomena do not contain such ex- 
treme pathological cases. 




FIG. 15: Attracting heteroclinic cycles may prevent 
convergence of ergodic averages. 

For dissipative systems, one often requires a that 
the set of Birkhoff-regular initial condition^^, for 
which Fourier averages of continuous observables 
converge, is of positive volume, i.e., that the system 
preserves a physical measur^^ supported on the at- 
tractor. In those cases, computing quotient maps 
using Fourier averages will identify basins of attrac- 
tions of the attractors, instead of attractors them- 
selves, since the points in the basin will have the 
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same averages as the points on the attractor. 

A weh-known exampld^ of a dissipative system 
for which trajectory averages do not converge for al- 
most ah open sets of initial conditions contains an 
attractive heteroclinic cycle (see Fig. 15). Ergodic 



measures are (5-measures supported on each of the 
fixed points, consequently, any invariant measure 
will be singular with respect to Lebesgue. Trajec- 
tories that approach heteroclinics spend longer and 
longer times along each of the exterior fixed points, 
possibly in such a manner that the finite time aver- 
ages do not converge. Notice that, while the condi- 



tions of Theorem 27 are formally satisfied, they are 
almost vacuous for practical purposes, as only the 
subsets of heteroclinic orbits have well-defined aver- 
ages, yet the convergence set is still of full measure, 
due to singular nature of the measure conserved. 

From the authors' experience in applying the av- 
eraging technique to dynamical systems that model 
physical phenomena for practical analysis, finding 
an appropriate convergence set E is rarely a prob- 
lem. A more practically significant issue comes from 
the different rates of convergence of finite to infinite 
averages for points within S. 

The rate of convergence of finite limits 

N-l 

An fix) :=^E/°^"(^)' 

n=0 



and their Fourier counterparts (67) on the set S is 
not uniform. It is a classical result that the rates 
of convergence over periods of length N can range 
from exponentially fast, i.e., 0{e~^^) for A > 0, for 
trajectories approaching exponentially stable fixed 
points, to algebraic convergence A^~", where a > 
is arbitrarily small, see, e.g., Petersen, §3. 
From a practical perspective, in a lot of cases the 
situation does not look so bleak. Trajectories on pe- 
riodic orbits and in strongly mixing regions achieve 
the rates of convergence of 0{N~^) and 0(A^~^/^), 
respectively. The slopes < a < 1/2 are to be ex- 
pected near homoclinic and heteroclinic orbits, es- 
pecially if such orbits are embedded within zones 
of intermittency, where trajectories get entrained 
around marginally stable fixed points for long times, 
before eventually moving away from them.^^ Such 
zones appear, for example, in perturbed hamil- 
tonian and volume-preserving systems. Studies 
of volume-preserving systems, in theoretical^^ and 
computational^ contexts, have shown that such re- 
gions are small in area. 

To illustrate, we plotted convergence errors for the 
Chirikov Standard Map (61), simulated for e = 0.18 



where the state space contains both regular and 
mixing regions. This b rief analysis is similar to 
those in literature .'^'^S'^ As an indicator of conver- 
gence speed, we compared the ergodic average, i.e.. 



Fourier average with = 0, after A^^o = 10^ it- 
erates with the average after (significantly) shorter 
number N <C N^q of iterates. To avoid making con- 
clusions based on a single observable, a truncated 
set of Fourier harmonics 

fk{x,p) = ^ exp [i27T{kxX + kpp)] , 

ZTT 

for k = {kx^ kp) G [—5, 5]^ C Z^, was used; this is the 
type of set used to practically approximate the quo- 



tient maps described in Section IV B and for each 
trajectory the largest absolute difference in averages 
over that set was taken as an indication of the con- 
vergence error. 



Figure 16a shows the error in averag es after a fixed 
time N = 10^ iterates, while Figures 16b, 16c, and 



16d show the time required for the error to stay 
within 10~^ for 100 iterates. It is clear that the 
speed of convergence varies across the state space, 
with convergence in chaotic regions slower than in 
regular regions. We proposed^ that the simulated 
time is varied for individual trajectory based on the 
relative convergence error, which resulted in effi- 
cient simulation runs with only a small subset of ini- 
tial conditions, in regions of intermittency, requiring 
long simulation times. 

This section dealt with the problem of extending 
finite-time averages into the infinite limit, to be able 
to approximate the limiting measures of empirical 
distributions. It might be surprising that even ex- 
plicitly finite time averages find their use in practical 
applications. In the next section, we demonstrate 
how to formulate continuous indicators of mixing 
and ergodicity, and use them in design of feedback 
control for technical systems. 



V. CONTINUOUS INDICATORS OF ERGODICITY 
AND MIXING 

When dynamical systems are analyzed as 
measure-preserving transformations, ergodicity and 
mixing property are among the first concepts dis- 
cussed in introductory textbooks .f^^I^ Let measure 
/i be preserved by the system. In plain language, er- 
godicity with respect to /i means that sets left invari- 
ant by the transformation/fiow are either full mea- 
sure /i, or /i- negligible. The mixing property, which 
implies ergodicity, means that any set of positive /i- 
measure will be distributed by the fiow according to 
the mixing measure /i. 

While the standard introductions to these 
properties'^ often include equivalence theorems that 
provide alternative formulations of ergodicity and 
mixing, cf. (60) and (71), all the definitions treat 



them as binary indicators: either the system is er- 
godic/mixing, or it is not. From the perspective of 
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(a) Difference in averages after N = 10^ and Noo — 10^ 
iterates (log^Q-scale). The bottom of the color scale 
includes all the values below —6. 
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(c) Number of steps required for the difference in 
averages compared to 10^ iterates to reach 10""^. 
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(b) Distribution of number of initial conditions over 
iterations required for the difference in averages 
compared to A^oo = 10^ iterates to reach 10~^. 
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(d) Number of steps required for the difference in 
averages compared to 10^ iterates to reach 10~^ 
(enlarged). 



FIG. 16: Convergence of ergodic averages. The difference in averages of Fourier harmonics for trajectory 
starting at x G M is computed as max^ \A^j}z{x) — A^^f]^[x)\ over /c G [— 5, 5]^ C . Step A^oo = 10^ 
was taken as "true" infinity. A total of 2025 initial conditions, seeded from a uniform rectangular grid, was 
used. Colors were interpolated from the first 500 iterates of each trajectory. 



applied dynamical systems, especially in the context 
of design of dynamics, having a binary indicator of a 
desired property is insufficient. Designers prefer to 
work with continuous indicators, e.g., an ergodicity 
indicator that takes values in [0, 1] where the ex- 
tremum would imply classical ergodicity, and the 
higher values would indicate how far, in some sense, 
the system is from being ergodic. The continuous in- 
dicators of ergodicity and mixing can be constr ucted 
from averages of functions along trajectories! ^^ * ^^ * 
The averaging process corresponds to an invariant 
measure, which is then compared to the a priori 
measure ja to infer how close the system is to be- 
ing /i-ergodic or /i-mixing. 

Designing systems to be ergodic or mixing has a 
number of practical applications. Consider the case 
of micro-mixers, devices whose task is to mix two 
fluids, react ants, such that when the reaction is ini- 



tiated, there are no "pockets" or unused reactants, 
and the reaction occurs uniformly in the vessel. At 
usual macro-scales, the fluids are easily mixed just 
by shaking them; on micro- and nano- scales this 
is not possible. Instead, micro-mixers use dynam- 
ical, advective transport, to reach the mixed state. 
An example where ergodicity is desirable comes from 
search-and-rescue missions, where helicopters or air- 
planes are used to scan a large area for survivors of 
airplane crashes and capsized boats. It is not triv- 
ial to design a flight path that ensures that the en- 
tire area is covered and that particular zones are 
searched more often or more thoroughly: such task 
can be phrased as design of an ergodic flight path. 
Both of these problems can be addressed by casting 
them into a framework of dynamical systems and re- 
quiring that the trajectories are mixing or ergodic. 
While trajectory averages were treated as a 
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FIG. 17: For finite time averaging functionals 

can be represented by empirical probability 
measures, whose distributions Cx,n can be thought 
of as sequences of S distributions supported on 
trajectories. 



technique for computation of projections onto 
eigenspaces of the Koopman operator, which re- 
quired their extension into infinity, in this section 
we will show what can be extracted from finite-time 
averages ( [67[ ), with = 0, i.e., the finite-time ver- 
sion of ergodic averages. The average of a bounded 
observable / : M ^ C over a finite-time trajectory 
segment, f ^ jj^ Yln=o f ° T^{x), is a linear con- 
tinuous functional. By the Riesz representation the- 
orem, the finite- average functionals are represented 
by empirical measures^ whose distributions, formally 
defined as 



with a test function / yields 

N-1 

= mT. S\p-T"{x)]f{p)dp 

N-1 

If, as an observable, we take the characteristic func- 
tion of a measurable set A, rx^N{A) = {cx,n^Xa) 
is the residence time: the fraction of time [0, N] 
that trajectory T^{x) spends inside set A. There- 
fore, the empirical measures {cx^n^ ') induced by av- 
eraging functionals are probability measures, for any 
initial condition x G M and N > {). For any finite 
time, empirical distributions are not absolutely con- 
tinuous with respect to the Lebesgue measure. In 
the limit N ^ oo, however, they converge in weak-* 
topology to ergodic measures;^^ their limits might be 
point masses, on equilibria, singular measures, e.g., 
on periodic orbits, or absolutely continuous mea- 
sures, e.g., in mixing regions. 

Ergodicity of the flow with respect to a flow- 
invariant measure fi can be stated as the condition 
that for [/i]-a.e. x G M, 

lim {cx,n,Xa) = /i(^), (71) 



N-1 



(70) 



can be thought of as strings (or ribbons) of ^- 
distributions, supported along the orbits (see Fig. 
17). 



The consequence of defining the empirical distri- 
bution using the Riesz representation theorem is the 
well-known equality between spatial and time aver- 
ages. With a bit of breadth in notation, the ex- 
pression (/, ^) = f{p)g{p)dp can be used to cou- 
ple functions / and measures g{p)dp. In the case 
of measures with density functions, i.e., absolutely 
continuous measures, both / and g can be taken 
as functions; when g{p)dp = d/jj{p) corresponds to a 
more general measure /i, / is taken as a test function, 
and g is in the class of distributions, or generalized 
functions, e.g., for a point mass /i, ^ is a Dirichlet 
6 distribution. Coupling the empirical measure ( 70 ) 



for all measurable sets A. The state-space points 
X e M whose trajectories satisfy the ergodicity con- 
dition are called /i-generic. The first step towards a 
continuous indicator of ergodicity was presented by 
Scott et al.!^, who used a Haar wavelet basis as ob- 
servables: these wavelets directly relate to indicator 
functions of a basis for open sets on a rectangular 
state space. Since wavelet bases are parametrized 
by scale at which sets are resolved, the constructed 
ergodic indicator, termed ergodicity defect^ acted as 
a coarse-grained version of ergodicity, reflecting an 
engineering perspective: instruments have finite pre- 
cisions, therefore, if the equalities (71) are satisfied 



for all sets coarser than precision of the instrument, 
the system can be thought of as ergodic for the par- 
ticular application. 

Based on a similar idea, Mathew and Mezic con- 
structed a different continuous ergodicity indicator. 
For declaring ergodicity, it is sufficient to verify the 
condition (71) on a basis for the Borel algebra, e.g., 
euclidean balls B{x,r) = {p G M : — x||2 < r}, 
whose indicator functions we label X{p,r)- Integrat- 
ing the deviations between trajectory averages and 
measures of sets, we obtain the empirical ergodicity 
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r)]\ dpdr, 
(72) 



where R is selected such that the largest ball in- 
cludes the entire state space M (cf. pseudodistance 
(|64|). It is almost immediate that the ergodicity is 
equivalent to lim^^oo Ex (n) = if x is selected as 
a /i- generic point. Strictly speaking, E^ compares 
the N ^ oo limit of empirical measures to the prior 
/i; to assure that the same behavior occurs almost 
everywhere, E^ can be integrated along the state 
space. 

The practical value of Ex{n) is in how its con- 
stituents are computed: {cx,n^Xip,r)) is computed 
as a time- average of an indicator function, which is 
very easily computed as the system is being simu- 
lated from an ODE. When /i is the volume measure, 
quantities /i[5(x,r)] are just volumes of euclidean 
balls, and are computed using well-known formulas. 
Consequently, evaluation of Ex{n) is simple, as it 
requires only evaluation of finite-time averages. 

On the other hand, fi can be a more detailed 
measure, as in the mentioned case of probabilities 
of target detect ion. ^^^I^ In those cases, evaluating 
(72) might be a more complicated effort, requiring 



careful spatial gridding to control errors in the in- 
tegral. Instead of measures of spherical sets, the 
metric Ex{n) can be expressed using Fourier coef- 
ficients of /i, facilitating the control of spatial scale 
resolution. The integrand in ( 72 ) can be interpreted 



as a difference between generalized expansion coeffi- 
cients, where the basis set is the set of characteristic 
functions, whose supports are a base for Borel sets: 

fi[B{p,r)] = (5/i,X(p,r)) . 

The djii stands for the (formal) density of the prior 
measure /i. 

If the basis X(p,r) is replaced by a harmonic ba- 
sis //e, then := {dja^fk) are just spatial Fourier 
coefficients, easily computed to very high orders by 
Fast Fourier Transform if we know the density d/n ex- 
plicitly. At the same time, practicality of evaluation 
of {cx,n^fk) is not sacrificed. The resulting empir- 
ical ergodicity is given by a metric induced by the 
negative-index Sobolev norm ||.||2 _5 on the space of 
distributions (cf. ([62|)), if both d/i and Cx^n 

are in it. 



\\Cx,n - d/J.\\ 



2,-s 



2^ n ^ fo^ whw \2]s 



kez- 



[l + i2n\\k\\,y 



(73) 



where the state space is, for simplicity, M T^, 
k G Z^. Wavevectors k G parametrize harmonic 
functions fk{x) = (27r)^/2gi27r/c.x^ 

When the order of the Sobolev space s is chosen 
as 5 = (D l)/2, \\cx,n Ex{n) are 

equivalentP^ there exist constants a > and f3 > 
such that 

<^ \\Cx,n - d/J^\\2_s < Ex{n) < p \\Cx^n " ^M||2 -s ' 

at all n. As a consequence, decay oi \\cx^n — dfiW^ _^ 
can be used as the proxy for computing decay of 
Exin) to detect ergodicity, since \\cx^n — djiW^ _g is 
easier to numerically evaluate for most measures fi. 

A practical application of Ex{n) can be seen on 
a model search-and-rescue problem, where an Un- 
manned Aerial Vehicle (UAV), e.g., a small he- 
licopter, would search an area containing a tar- 
get whose position is estimated by a probability 
densit}'^^. The search path for the vehicle can be 
planned using a dynamical system: the searcher is 
treated as a passive particle, possibly with a non- 
zero inertia, in a fiuid-like fiow. The probability 
density is a product of gaussians, modeling a priori 
estimates of the target, and discontinuous indica- 
tor functions, modeling foliage on the ground where 
the UAV has no visibility. The goal is to design 
the dynamics of the fiow, such that the trajectory 
traced out by a particle avoids foliage and explores 
the feasible area according to probability of finding 
the target. 

The devised control algorithm used Ex{t) (in 
continuous-time setting) as the optimization func- 
tion in a Hamilton-Jacobi-Bellman framework. In- 
stead of optimizing for decay of Ex{t) over a fi- 
nite horizon t G [0,T], a greedy approach is chosen, 
where the optimization horizon is shrunk to the in- 
stance by taking T ^ 0, resulting in a closed- form 
expression feedback law which relies on finite-time 
averages of harmonic observables, computed along 
the path that the searcher traveled. As a result, 
trajectories designed for a team of UAVs searched 
the area efficiently, with minimal crossings over the 
zones with high foliage (Fig.[l8|, with further details 
presented in^^ 

In addition to quantifying ergodicity, the W^'~^ 
norm can be used to quantify mixing with respect 
to a target density in which context it was 
originally developed. In the final design of a 
micro- mixer, an initial concentration of fiuid reac- 
tants p was advected by a dynamical system, whose 
effect was captured by the Perron-Frobenius opera- 
tor P^p. To quantify mixing, the metric M(n) := 
||P^p — dp\\2 _g was computed with the order s = 
1/2. By optimizing the dynamics of the ffuid ffow, a 
quick mixing of react ants was successfully achieved 
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(a) Target probability 
distribution and initial searcher 
positions. Probability density is 
positive and constant on white 
region, zero on gray. 



(b) Paths of searchers, showing 
the searchers sampling the area 
with positive density. 



(c) Decay of | 



■ d^W^ _s (labeled by $ in 



the original paper). 



FIG. 18: Searching for a target by optimizing the decay of empirical ergodicity Ex{r^ through 
computation of \\cx^n — djiW^ _g for several trajectories. (Original in Mathew and Mezic^, Physica D: 
Nonlinear phenomena by North-Holland. Reproduced with permission of North-Holland in the format 
reuse in a journal/magazine via Copyright Clearance Center.) 



(Fig. 19), both in numerical and experimental set- 
tings. 

The main difference between M(n) and Ex{n) is 
that the mixing metric compares instantaneous ad- 
vected density P^p to the prior /i, which, in the dual 
Koopman framework, is equivalent to using U'^f. 
Recall that the ergodicity metric used the temporal 
averages ^ T.^~o[U''f] (cf. ([TO])). The other differ- 
ence between mixing and ergodicity indicators can 
be interpreted as the requirement on the scales at 
which measure /i is sampled: the ergodicity indica- 
tor Ex{n) gives more weight to larger spherical sets, 
while in M(n), differences over all sets are weighted 
the same, regardless of the set size, due to rescaling 
of the measure dxdr by the measure of the spheri- 
cal set to dxdr/ii{x^ r). Therefore, those trajectories 
that lead to uniform decay of Ex{n) will first appear 
well distributed in /i over larger sets, and only then 
will they sample /i on smaller scales. Conversely, 
M{n) does not make such a distinction, requiring 
that the trajectories distribute on all scales equally 
fast. 

In practice, the choice between Ex{n) and M{n) 
is driven by the application: for mixing of fluids, it 
was important that two reactants were in contact at 
the instant of reaction initiation, e.g., initiation of 
burning in the chamber. In the search-and-rescue 
application, instantaneous behavior was not of the 
essence, as the goal was to find the target over the 
entire course of the search mission, not necessarily 



to have the same (high) probability to locate the 
target at any particular time instance. In future 
applications, a similar analysis of the problem would 
decide which property, ergodicity Ex{n) or mixing 
M{n), would be a more appropriate design criterion. 

The indicators for ergodicity and mixing presented 
in this section are based on representations of er- 
godic measures in the space of Fourier coefficients. 
They therefore represent ergodic measures as accu- 
rately as it is possible through their Fourier coef- 
ficients: numerically, approximations will converge 
quicker the smoother the ergodic measure studied. 
However, it is well known that the shape of the er- 
godic measure's density does not determine uniquely 
the dynamical system. In this sense, the indicators 
are not intrinsic to the dynamical system studied. 
For example, if one would use a unique ergodic mea- 
sure /i for a map Ti and seek to drive the ergodic 
measure of another system T2 to match /i, there 
would be no guarantee that trajectories of Ti and 
T2 would be conjugate, without further restrictions 
on properties of maps Ti and T2. 



VI. CONCLUSIONS 

We have reviewed the theoretical aspects and ap- 
plications of the spectral theory of the Koopman op- 
erator in dynamical systems. The use of these con- 
cepts holds promise to provide a theory that extends 
and complements tools from geometrical dynamical 
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FIG. 19: Simulation of mixing of react ants in a micro-chamber using dynamics designed by optimizing a 
W~^'^ distance between advected density and uniform target density. (Original in Mathew et al.l^, 
Journal of fluid mechanics by Cambridge University Press. Reproduced with permission of Cambridge 
University Press in the format reprint in a journal via Copyright Clearance Center.) 



systems theory that enabled so much development 
in science and technology over the last century. The 
presented material described three parallel branches 
of Koopman operator analysis: (i) Koopman mode 
analysis, (ii) eigenquotient analysis, and (iii) indica- 
tors of ergodicity and mixing. The Koopman modes 
generalize the notion of linear eigenmodes, known 
from, e.g., linear mechanical vibration theory, to 
the nonlinear context, without linearizing the dy- 
namics first. It is interesting that projections onto 
eigenspaces of the Koopman operator are achieved 
via an extension of Laplace and Fourier analysis, and 
that such extension works in the nonlinear case. The 
Koopman mode analysis has proved especially use- 
ful in analysis of dynamics in infinite-dimensional 
state spaces, which were observed using a high- 
dimensional measurements, e.g., thermal dynamics 
of a building system, with a distributed measure- 
ment of temperatures or fluids experiments. 

The eigenquotient analysis generalizes the analy- 
sis of smooth integrals of motion, which is central 
in theoretical mechanics, to systems which do not 
have any smooth invariants. The non-smooth eigen- 
functions are constructed using infinite-time averag- 
ing along trajectories of the dynamical system. The 
eigenfunctions are used to construct a geometry for 
families of invariant sets which can be used to ex- 
tract invariant regions that locally resemble phase 
portraits of integrable Hamiltonian systems. The 
same analysis can not only be applied to the study 
of invariant sets, but also periodic and wandering 
sets. 

Finally, we showed how the trajectory averages 
can be used even when only finite-time data is avail- 
able. The values of finite- time averages of functions 
along trajectories can be used to evaluate quantita- 



tive indicators of how closely the system is to ac- 
curately sampling a prior measure. Such indicators 
can be used as optimization criteria in a feedback al- 
gorithm for trajectory planning of Unmanned Aerial 
Vehicles. 

The theory reviewed here was introduced and de- 
veloped mostly for measure-preserving determinis- 
tic systems, including dynamics on the attractor 
of a dissipative system. Much more work remains 
to be done in the case of dissipative systems and 
extensions to non-smooth and hybrid (determinis- 
tic/stochastic) case. 
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